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
Radial Release
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] w0 = 5; % [rad/s] 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 = Ig_inv*Hgx % [rad/s]
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
Small Satellites 