Universidad Europea de Madrid

School of Architecture, Engineering and

Design

Mid-Term Problem Set

Hans Lluis Herrmann Jane

Department of Aerospace Engineering

Course: Space Vehicles and Missiles

November 21, 2023

 Exercise 1
In this exercise, I present my approach and solution to the problem of deter-
mining the position and velocity of the International Space Station (ISS) after
72 hours, given its initial position and velocity. This exercise demonstrates the
application of orbital mechanics principles using MATLAB.
The motion of an object in orbit around the Earth is described by the two-body
problem. This simplification allows us to use Newton’s law of universal gravi-
tation and his second law of motion to derive the equations of motion for the ISS.

I used MATLAB to numerically integrate the equations of motion over a pe-

riod of 72 hours. This approach is based on the assumption that the ISS follows
a Keplerian orbit, which is a good approximation for short-term predictions.
The following MATLAB code was used to compute the new position and velocity
vectors:
1 % Given data
2 R0 = [3781.9; 2444.1; 5119.1]; % Initial position vector ( km )
3 V0 = [ -2.2634; 7.0792; -1.7514]; % Initial velocity vector ( km / s )
4 t0 = 0; % Initial time ( hours )
5 t_end = 72; % Time after which to determine the state ( hours )
6
7 % Constants
8 mu_earth = 398600; % Earth gravitational parameter ( km ^3/ s ^2)
9
10 % Runge - Kutta integration
11 options = odeset ( ’ RelTol ’ , 1e -8 , ’ AbsTol ’ , 1e -8) ; % Set tolerances
12 [t , state ] = ode45 ( @ (t , y ) o r b i t a l _ d y n a m i c s (t , y , mu_earth ) , [ t0 ,
t_end ] , [ R0 ; V0 ] , options ) ;
13
14 % Extracting position and velocity vectors after 72 hours
15 R_final = state ( end , 1:3) ’;
16 V_final = state ( end , 4:6) ’;
17
18 disp ( ’ Position vector after 72 hours : ’) ;
19 disp ( R_final ) ;
20 disp ( ’ Velocity vector after 72 hours : ’) ;
21 disp ( V_final ) ;
22
23 % Plotting the orbits
24 figure ;
25 plot3 ( state (: , 1) , state (: , 2) , state (: , 3) , ’ LineWidth ’ , 2) ;
26 grid on ;
27 axis equal ;
28 title ( ’ Orbit of the International Space Station ’) ;
29 xlabel ( ’X ( km ) ’) ;
30 ylabel ( ’Y ( km ) ’) ;
31 zlabel ( ’Z ( km ) ’) ;
32
33 % Function for orbital dynamics
34 function dydt = o r b i t a l _ d y n am i c s (~ , y , mu )
35 % y = [ Rx ; Ry ; Rz ; Vx ; Vy ; Vz ]
36
37 % Extracting position and velocity vectors

1
38 R = y (1:3) ;
39 V = y (4:6) ;
40
41 % Calculating acceleration
42 r = norm ( R ) ;
43 acceleration = - mu / r ^3 * R ;
44
45 % Forming the state derivative
46 dydt = [ V ; acceleration ];
47 end

The final position and velocity of the ISS were calculated and visualized in
a 3D plot, showing the trajectory over 72 hours.
The results provide a reasonable approximation of the ISS’s position and velocity
after 72 hours, given the initial conditions and the assumptions made in the
model.

2
 Exercise 2
This exercise involves determining the position of the International Space Sta-
tion (ISS) 46.493 minutes after the epoch time given in a Two-Line Element
(TLE) set. The TLE provides essential orbital parameters that are used to
compute the orbit of the ISS.

A TLE is a data format encoding a list of orbital elements of an Earth-

orbiting object for a given point in time. The ISS’s orbit can be determined
using these elements, which include the inclination, right ascension of the as-
cending node (RAAN), eccentricity, argument of perigee, mean anomaly, and
mean motion. The orbit is propagated forward in time to find the position of
the ISS after a specified duration.

The MATLAB code reads the TLE data, extracts the necessary orbital ele-
ments, and then uses these elements to compute the ISS’s position after 46.493
minutes. The calculation involves solving Kepler’s equation for the eccentric
anomaly and then determining the true anomaly. The orbital position is trans-
formed from the orbital plane into a 3D coordinate system.
1 % ISS TLE Data
2 TLE1 = ’1 25544 U 98067 A 21327. 09520073 -.00003894 00000+0 -63390 -4
0 9991 ’;
3 TLE2 = ’2 25544 51.6433 276.9828 0005100 229.5055 237.2561
1 5 . 4 8 6 1 3 1 3 2 3 1 3 2 2 1 ’;
4
5 % Extracting orbital elements from TLE
6 inclination = str2double ( TLE2 (9:16) ) ; % Inclination [ degrees ]
7 RAAN = str2double ( TLE2 (18:25) ) ; % Right Ascension of the Ascending
Node ( RAAN ) [ degrees ]
8 eccentricity = str2double ([ ’ 0. ’ , TLE2 (27:33) ]) ; % Eccentricity
9 a r g u m e n t O f P e r i g e e = str2double ( TLE2 (35:42) ) ; % Argument of Perigee
[ degrees ]
10 meanAnomaly = str2double ( TLE2 (44:51) ) ; % Mean Anomaly [ degrees ]
11 meanMotion = str2double ( TLE2 (53:63) ) ; % Mean Motion [ revolutions
per day ]
12
13 % Convert mean motion to semi - major axis
14 mu = 398600; % Standard gravitational parameter for Earth [ km ^3/ s
^2]
15 r e v s P e r D a y T o R a d P e r S e c = 2 * pi / (24 * 3600) ; % Convert revolutions
per day to radians per second
16 n = meanMotion * r e v s P e r D a y T o R a d P e r S e c ; % Mean motion [ rad / s ]
17 a = ( mu / n ^2) ^(1/3) ; % Semi - major axis [ km ]
18
19 % Calculate the future position and velocity (46.493 minutes later )
20 deltaTime = 46.493 * 60; % Convert minutes to seconds
21 [ futurePosition , futur eVeloci ty ] = p ropagate Orbit (a , eccentricity ,
inclination , RAAN , argumentOfPerigee , meanAnomaly , deltaTime ,
mu ) ;
22
23 % Plot the future position of the ISS
24 figure ;

3
25 plot3 ( futur ePositio n (1) , fu turePosi tion (2) , futureP osition (3) , ’o ’)
;
26 title ( ’ Future Position of the ISS ’) ;
27 xlabel ( ’X ( km ) ’) ;
28 ylabel ( ’Y ( km ) ’) ;
29 zlabel ( ’Z ( km ) ’) ;
30 grid on ;
31
32 % Function to propagate the orbit
33 function [ position , velocity ] = propag ateOrbit (a , e , i , Omega ,
omega , M0 , deltaTime , mu )
34 % Convert to radians
35 i = deg2rad ( i ) ;
36 Omega = deg2rad ( Omega ) ;
37 omega = deg2rad ( omega ) ;
38 M0 = deg2rad ( M0 ) ;
39
40 % Solve the eccentric anomaly
41 M = M0 + sqrt ( mu / a ^3) * deltaTime ;
42 E = k e p l e r E q u a t i o n S o l v e r (M , e ) ;
43
44 % Calculate position and velocity in perifocal coordinates
45 r = a * (1 - e * cos ( E ) ) ; % Radial distance
46 x = r * cos ( E ) - a * e ; % X position
47 y = r * sin ( E ) * sqrt (1 - e ^2) ; % Y position
48
49 % Velocity in perifocal coordinates
50 vx = - sqrt ( mu / a ) * sin ( E ) ;
51 vy = sqrt ( mu / a ) * ( cos ( E ) - e ) * sqrt (1 - e ^2) ;
52
53 % Rotate to inertial coordinates
54 position = rotateCoords ([ x ; y ; 0] , Omega , i , omega ) ;
55 velocity = rotateCoords ([ vx ; vy ; 0] , Omega , i , omega ) ;
56 end
57
58 % Function to solve Kepler ’ s equation
59 function E = k e p l e r E q u a t i o n S o l v e r (M , e )
60 % Iterative method to solve Kepler ’ s equation
61 E = M;
62 for j = 1:10
63 E = E - ( E - e * sin ( E ) - M ) / (1 - e * cos ( E ) ) ;
64 end
65 end
66
67 % Function to rotate coordinates
68 function rotatedCoords = rotateCoords ( coords , Omega , i , omega )
69 % Rotation matrix
70 R3_Omega = [ cos ( Omega ) , sin ( Omega ) , 0; - sin ( Omega ) , cos ( Omega ) ,
0; 0 , 0 , 1];
71 R1_i = [1 , 0 , 0; 0 , cos ( i ) , sin ( i ) ; 0 , - sin ( i ) , cos ( i ) ];
72 R3_omega = [ cos ( omega ) , sin ( omega ) , 0; - sin ( omega ) , cos ( omega ) ,
0; 0 , 0 , 1];
73
74 % Rotate the coordinates
75 rotatedCoords = R3_Omega ’ * R1_i ’ * R3_omega ’ * coords ;
76 end

4
 The MATLAB code successfully parses the TLE data and computes the ISS’s
position in the Earth-Centered Inertial (ECI) coordinate system at the target
time. The resulting position vector provides a snapshot of the ISS’s location in
its orbit after the given time advancement. Here is the representation of the
results in a visualized way:

5
 Exercise 3
In this exercise, the task is to determine the position and velocity of a spacecraft
orbiting the Moon, given its orbital elements, after a duration of 27.3217 days.
The problem incorporates the effects of the Moon’s J2 perturbation.

The spacecraft’s orbit is characterized by its semimajor axis, eccentricity, in-

clination, right ascension of the ascending node (RAAN), argument of perigee,
and true anomaly. The gravitational influence of the Moon, including its J2
perturbation (which accounts for the oblateness of the Moon), affects the orbit.
The J2 term causes the orbital parameters to change over time, particularly the
RAAN and argument of perigee.

The MATLAB code computes the initial position and velocity of the space-
craft from the given orbital elements. It then propagates the orbit over 27.3217
days, taking into account the Moon’s J2 perturbation. This requires solving
the orbital dynamics equations using numerical integration.

1 % Define the given orbital elements

2 a = 10500; % Semimajor axis in km
3 e = 0.1; % Eccentricity
4 i = 55 * pi /180; % Inclination in radians
5 RAAN = 238.4 * pi /180; % Right Ascension of Ascending Node in
radians
6 omega = 31.6 * pi /180; % Argument of perigee in radians
7 theta = 65.4 * pi /180; % True anomaly in radians
8
9 % Gravitational parameter of the Moon
10 mu = 4903; % in km ^3/ s ^2
11
12 % J2 perturbation term
13 J2 = 202.7 e -6;
14
15 % Equatorial radius of the Moon
16 R_moon = 1738.1; % in km
17
18 % Convert orbital elements to position and velocity vectors
19 [p , v ] = orb2rv (a , e , i , RAAN , omega , theta , mu ) ;
20
21 % Time to propagate
22 T = 27.3217 * 24 * 3600; % Time in seconds
23
24 % Propagate the orbit using ode45 with J2 perturbation
25 options = odeset ( ’ RelTol ’ , 1e -8 , ’ AbsTol ’ , 1e -8) ;
26 [t , y ] = ode45 ( @ (t , y ) twobodyJ2 (t , y , mu , R_moon , J2 ) , [0 T ] , [ p ;
v ] , options ) ;
27
28 % Extract the final position and velocity
29 final _positi on = y ( end , 1:3) ;
30 final _veloci ty = y ( end , 4:6) ;
31
32 % Display results
33 disp ( ’ Final position vector ( km ) : ’) ;

6
34 disp ( fin al_posi tion ) ;
35 disp ( ’ Final velocity vector ( km / s ) : ’) ;
36 disp ( fin al_velo city ) ;
37
38 % Plotting the orbit
39 figure ;
40 plot3 ( y (: ,1) , y (: ,2) , y (: ,3) , ’b ’) ;
41 hold on ;
42 plot3 ( final _positio n (1) , fi nal_posi tion (2) , final_p osition (3) , ’ ro ’
);
43 title ( ’ Orbit of the Spacecraft around the Moon ’) ;
44 xlabel ( ’x ( km ) ’) ;
45 ylabel ( ’y ( km ) ’) ;
46 zlabel ( ’z ( km ) ’) ;
47 grid on ;
48
49 % Supporting functions
50 function [p , v ] = orb2rv (a , e , i , RAAN , omega , theta , mu )
51 % This function converts orbital elements to position and
velocity vectors
52
53 % Perifocal coordinates
54 p_pf = ( a * (1 - e ^2) ) / (1 + e * cos ( theta ) ) ;
55 r_pf = [ p_pf * cos ( theta ) ; p_pf * sin ( theta ) ; 0];
56 v_pf = sqrt ( mu / ( a * (1 - e ^2) ) ) * [ - sin ( theta ) ; e + cos ( theta
) ; 0];
57
58 % Rotation matrices
59 R3_W = [ cos ( - RAAN ) sin ( - RAAN ) 0; - sin ( - RAAN ) cos ( - RAAN ) 0; 0 0
1];
60 R1_i = [1 0 0; 0 cos ( - i ) sin ( - i ) ; 0 - sin ( - i ) cos ( - i ) ];
61 R3_w = [ cos ( - omega ) sin ( - omega ) 0; - sin ( - omega ) cos ( - omega ) 0;
0 0 1];
62
63 % Transform to inertial frame
64 p = R3_W * R1_i * R3_w * r_pf ;
65 v = R3_W * R1_i * R3_w * v_pf ;
66 end
67
68 function dy = twobodyJ2 (~ , y , mu , R , J2 )
69 % This function calculates the acceleration due to gravity and
J2 perturbation
70
71 % Position and velocity
72 r = y (1:3) ;
73 v = y (4:6) ;
74
75 % Distance from the Moon
76 r_norm = norm ( r ) ;
77
78 % Acceleration due to gravity
79 a_grav = - mu / r_norm ^3 * r ;
80
81 % J2 perturbation
82 x = r (1) ;
83 y = r (2) ;
84 z = r (3) ;

7
85 a_J2 = -(3/2) * J2 * ( mu * R ^2 / r_norm ^5) * [ x * (5* z ^2/ r_norm
^2 - 1) ; y * (5* z ^2/ r_norm ^2 - 1) ; z * (5* z ^2/ r_norm ^2 - 3) ];
86
87 % Total acceleration
88 a = a_grav + a_J2 ;
89
90 % Output derivative of state
91 dy = [ v ; a ];
92 end

The code outputs the position and velocity vectors of the spacecraft in the
Moon-centric frame after 27.3217 days. The results reflect the combined effect of
the Moon’s gravitational field and its oblateness, which is significant for accurate
long-term orbital predictions.

8
Exercise 4
Given two points on a geocentric orbit, we have the following data:
• Altitude z1 = 2555 km and true anomaly θ1 = 120◦ .
• Altitude z2 = 925 km and true anomaly θ2 = 60◦ .
The objective is to find the orbit’s eccentricity (e), the altitude of perigee (zp ),
the semimajor axis (a), and the orbital period (T).
The calculated orbital parameters are:
• Eccentricity e = 0.20096
• Altitude of Perigee zp = 314.57 km

• Semimajor Axis a = 8367.01 km

• Orbital Period T = 7616.70 seconds
The radial distance r from the center of the Earth to the satellite is given
by the altitude plus Earth’s radius. The equation for r in terms of the specific
angular momentum (h) and the eccentricity (e) is:
h2
µ
r= (1)
1 + e cos(θ)
where µ is the standard gravitational parameter and θ is the true anomaly.
Given r1 , r2 , θ1 , and θ2 , we can set up a system of equations to solve for e and
h.
The altitude of perigee zp is found by:

zp = rp − Rearth (2)
where rp is the perigee distance, calculated using the eccentricity:

h2 /µ
rp = (3)
1+e
The semimajor axis a can be found from:

h2
a= (4)
µ(1 − e2 )
Finally, the orbital period T is determined using Kepler’s third law:
s
a3
T = 2π (5)
µ

9
Exercise 5
This exercise involves calculating the specific orbital energy of two spacecraft in
different orbits around the Earth. The first is in a circular orbit, and the second
in an elliptical orbit.
The specific orbital energy (ϵ) in an orbit is given by:

v2 µ
ϵ= − (6)
2 r
where v is the orbital velocity, µ is the standard gravitational parameter of
Earth, and r is the distance from the Earth’s center to the spacecraft.
For a circular orbit, v is calculated using:
r
µ
v= (7)
r

For an elliptical orbit, the specific orbital energy is constant and given by:
µ
ϵ=− (8)
2a
where a is the semi-major axis of the orbit.
The specific energy for a spacecraft in a circular orbit with a radius of
10,000 km was calculated using the formula for circular orbits. Whereas,
the specific energy for a spacecraft in an elliptical orbit at perigee (altitude
of 5000 km) was calculated using the formula for elliptical orbits, taking
the average of the perigee and apogee distances as the semi-major axis.

- The specific orbital energy for the circular orbit is approximately −19, 930, 000 J/kg.
- The specific orbital energy for the elliptical orbit at perigee is approxi-
mately −12, 173, 966 J/kg.

10