Home » Space Flight/Orbital Mechanics Excercise

# Category Archives: Space Flight/Orbital Mechanics Excercise

## Bi-Elliptic Hohmann Transfer

In this example we calculate the total change in speed required for a bi-elliptic Hohmann transfer from a geocentric circular orbit of 7200 km radius to circular orbit of 125000 km radius. The apogee of the first transfer ellipse is 190000 km.

``` clc;
clear all;

R_i  = 7200;     % [km]
R1_a = 190000;   % [km]
R_f  = 125000;   % [km]

mu   = 398600;   % [km^3/s^2] Earth’s gravitational parameter
% For initial circular orbit
V_i = (mu/R_i)^0.5;
% Speed at apogee and perigee for the first transfer ellipse
V1_a = (2*mu*R_i/(R1_a*(R1_a+R_i)))^0.5;
V1_p = (2*mu*R1_a/(R_i*(R1_a+R_i)))^0.5;
% Semimajor axes of the first transfer ellipse
a1 = (R_i + R1_a)/2;
% Speed at apogee and perigee for the second transfer ellipse
V2_a = (2*mu*R_f/(R1_a*(R1_a+R_f)))^0.5;
V2_p = (2*mu*R1_a/(R_f*(R1_a + R_f)))^0.5;
% Semimajor axes of the second transfer ellipse
a2 = (R_f + R1_a)/2;
% For target circular orbit
V_f = (mu/R_f)^0.5;

% For bi-elliptic maneuver the total speed change required
dV =  abs(V_i - V1_p)+ abs(V1_a - V2_a) + abs(V_f - V2_p); %[km/s]

% Time required for transfer
t_bi = pi/(mu)^0.5*(a1^1.5+a2^1.5);     %[s]

fprintf('Total speed change = %4.4f [km/s]\n',dV);
fprintf('Time required for transfer = %4.2f [hours]\n',t_bi/3600);```
```Total speed change = 3.9626 [km/s]
Time required for transfer = 129.19 [hours]```

## Simple orbital mechanics problems

An unmanned satelite orbits the Earth

```clear all;
clc;
rp = 7000;                      %[km] Perigee radius
ra = 70000;                     %[km] Apogee radius
mu = 398600;                    %[km^3/s^2]

ecc = (ra-rp)/(ra+rp)           %Eccentricity of the orbit
a   = (ra+rp)/2                 %[km] Semimajor axis
T   = 2*pi/mu^0.5*a^1.5/3600    %[h] Period of the Orbit
E_spec = -mu/(2*a)              %[km^2/s^2]Specific energy```
```ecc =

0.8182

a =

38500

T =

20.8833

E_spec =

-5.1766```

Find true anomaly when Alt = 1000 km;

```R_earth = 6378;
theta = acos(a*(1 -ecc^2)/((R_earth+1000)*ecc)-1/ecc)*180/pi %[deg]```
```theta =

27.6069```

Velocity components at this anomaly

```hmu = (a*(1 -ecc^2)/mu)^0.5;
vr = 1/hmu *ecc*sin(theta)          %radial component
vt = 1/hmu*(1 + ecc*cos(theta))     %tangential component```
```vr =

2.8342

vt =

2.0001```

Velocities at Apogee and Perigee

```vp = 1/hmu*(1 + ecc)
va = 1/hmu*(1 - ecc)```
```vp =

10.1751

va =

1.0175```

The spacecraft is in a 250 km by 300km low eart orbit. How long does it take to coast from perigee to appogee. Semimajor axis of this eliptical orbit

```a = R_earth + (250+300)/2; % [km]
% Time which takes to fly from Perigee to Appogee equals to half of the
% Orbit period
Th = pi/mu^0.5*a^1.5/60 %[min]```
```Th =

45.0045```

The altitude of a satelite in an elliptical orbit around the earth is 1600 km at apogee and 600 km at perigee.

``` ra = R_earth + 1600;
rp = R_earth +600;
% Eccentricity of the Orbit
ecc = (ra - rp)/(ra + rp)
% Semimajor axis;
a = (ra + rp)/2
% The orbital Speed at perigee and apogee
hmu = (a*(1 -ecc^2)/mu)^0.5;
vp = 1/hmu*(1 + ecc)
va = 1/hmu*(1 - ecc)
% The period of the Orbit
T = 2*pi*a^1.5/mu^0.5/60 %[min]```
```ecc =

0.0669

a =

7478

vp =

7.8065

va =

6.8280

T =

107.2601```

A satellite is palced into an orbit at perigee at an altitude of 1270 km with a speed of 9 km/s. Calculate the flight path angle gamma and altitude of the satelite at a true anomaly of 100deg.

```altp  = 1270;    %[km]
vp    = 9;       %[km/s]
theta = 100;     %[deg]
rp = altp + R_earth;
h = rp * vp
ecc = h^2/(rp*mu) - 1
gamma = atan( ecc*sind(theta)/(1 + ecc*cosd(theta)));
gamma = 180/pi*gamma    %[deg]```
```h =

68832

ecc =

0.5542

gamma =

31.1256```

Altitude of the satelite at a true anomaly of 100 [deg].

```rd = h^2/mu*(1/(1 + ecc* cosd(theta)));
altd = rd - R_earth```
```altd =

6.7738e+003```

A satelite is launched into the orbirt at an altitude of 640 km with a speed 9.2 km/s and flight path angle of 10 deg. Calculate the true anomaly of the lunch poinnt and the period of the orbit

```alt = 640;           %[km]
rs  = R_earth + alt; %[km]
vel = 9.2;           %[km/s]
gamma = 10;          %[deg]
vp = vel*cosd(gamma);
vr = vel*sind(gamma);
h  = rs * vp;
ac = (h^2/(mu*rs)-1);
bc = (vr*h/mu);
theta = atan(bc/ac)*180/pi  %Anomaly of the lunch point
%Period of the orbit
%Eccentricity
ecc = ac/cosd(theta);
T = 2*pi/mu^2*(h/(1 - ecc^2)^0.5)^3  %[sec]```
```theta =

29.7830

T =

1.6075e+004```

A satellite has a perigee and apogee altitudes of 250 km and 42000km. Calculate the orbit period, eccentricity and the maximum speed.

```rp = R_earth + 250; %km
ra = R_earth + 42000;%km

%Semimajor axis of this eliptical orbit
a =( rp + ra )/2;
T = 2*pi/mu^0.5 * a ^(3/2)/60       %[min]
ecc = (ra - rp)/(ra+rp)
h = (mu*(1 + ecc)*rp)^0.5;
vp = h/rp                           %[km/s]```
```T =

756.5368

ecc =

0.7590

vp =

10.2852```

A rocket lunched from the surface of the earth has a speed 8.85 km/s when powered flight ends at an altitude of 550 km. The flight path angel at this time is 6 degree.Determinine eccentricity

```vel  = 8.85 ;           %[km/s]
ra = R_earth + 550;     %[km]
gamma = 6;              %[deg]
vd = vel * cosd(gamma);
h  = ra * vd;
ect = h^2/(mu*ra)-1;
vr = vel* sind(gamma);
est = h*vr/mu;
theta = atan(est/ect);
ecc = ect/cos(theta)```
```ecc =

0.3742```

Published with MATLAB® 7.10

## Euler rotation example, Rotation matrix, Quaternion, Euler Axis and Principal Angle

A classical Euler rotation involves first a rotation about e3 axis, then one about the e1 axis and finally a rotation about the e3 axis.

```% Compute the elements of rotation matrix R
clc;
psi   = pi/4;
theta = pi/4;
fi    = pi/6;

R3_psi   = [cos(psi) sin(psi) 0;
-sin(psi) cos(psi) 0;
0  0        1]
R1_theta = [1 0 0;
0 cos(theta)  sin(theta);
0 -sin(theta) cos(theta)]
R3_fi    = [cos(fi)   sin(fi) 0;
-sin(fi) cos(fi) 0;
0 0 1]```
```R3_psi =

0.7071    0.7071         0
-0.7071    0.7071         0
0         0    1.0000

R1_theta =

1.0000         0         0
0    0.7071    0.7071
0   -0.7071    0.7071

R3_fi =

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

Multiplying each rotation matrix

` R = R3_fi*R1_theta*R3_psi`
```R =

0.3624    0.8624    0.3536
-0.7866    0.0795    0.6124
0.5000   -0.5000    0.7071```

From Rotation Matrix Compute Quaternion

```q4 = 0.5*(1 + R(1,1)+ R(2,2) + R(3,3))^0.5;
q1 = (R(2,3) - R(3,2))/(4*q4);
q2 = (R(3,1) - R(1,3))/(4*q4);
q3 = (R(1,2) - R(2,1))/(4*q4);
q  = [q1 q2 q3 q4]
norm_q = norm(q)            % Checking that the norm of q = 1
qv     = [q1 q2 q3];```
```q =

0.3794    0.0500    0.5624    0.7330

norm_q =

1.0000```

Compute Euler axis and its components along the axis E1,E2,E3 unit vector. Compute the Euler principle angle.

```Euler_angle = 2*acos(q4)*180/pi            % [deg]
Euler_axis  = qv/sind(Euler_angle/2)       % [E1,E2,E3]
norm_E      = norm(Euler_axis)```
```Euler_angle =

85.7293

Euler_axis =

0.5577    0.0734    0.8268

norm_E =

1.0000```

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

A satellite is to be completely despin using a two-mass yo-yo device with tangential release. Assume the spin axis of moment of inertia of the satellite is C = 200 kg · m^2 and the initial spin rate w0 = 5 [rad/s]. The total yo-yo mass is 4 kg, and the radius of the spacecraft is 1 meter.

```clc;
C  = 200;       % [kg*m^2]
m  = 4;         % [kg]
R  = 1;         % [m]
K = 1 + C/(m*R^2); %     Nondimensional factor K```

The required cord length l

```l = R*(K)^0.5       %[m]
% Time t to despin
t = (K)^0.5/w0      %[s]
% Tension in the yo-yo cables
t = 0:0.001:5;
N = 2*K*C*w0^3*t./(R*(K + w0^2*t.^2).^2);
plot(t,N);
xlabel('time [s]');
ylabel('Tension N [N]');
grid on;
[maxN,i] = max(N);```
```l =

7.1414

t =

1.4283```

The maximum tension in the yo-yo cables [N];

``` maxN
% which ocurs at
t = i*0.001  % [s]```
```maxN =

454.7542

t =

0.8260```

The speed of the masses at release;

`v = R*w0*(K)^0.5  % [m/s]`
```v =

35.7071```

The angle rotated by the satellite during the despin;

`theta = (pi/2 - 1)*(K)^0.5   % [rad]`
```theta =

4.0763```

Cord length required for complete despin

`l = R*(K^0.5 - 1)            % [m]`
```l =

6.1414```

Published with MATLAB® 7.10

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