MAE 8 - Fall 2022 Project
Simulation of Satellites in Low Earth Orbit
1. Introduction
In this project, you are to simulate the trajectories of satellites in low Earth orbits. The
forces that act on the satellites are gravitational force and drag. From Newton’s second law,
the orbital motion of the satellites can be described by the following differential equations:
∂U
∂t
= −GMe X
(X2 + Y 2 + Z2)3/2
− CdρaAs
2m
U
√
U2 + V 2 +W 2,
∂V
∂t
= −GMe Y
(X2 + Y 2 + Z2)3/2
− CdρaAs
2m
V
√
U2 + V 2 +W 2,
∂W
∂t
= −GMe Z
(X2 + Y 2 + Z2)3/2
− CdρaAs
2m
W
√
U2 + V 2 +W 2, (1)
∂X
∂t
= U,
∂Y
∂t
= V,
∂Z
∂t
= W,
where t is time (in seconds); X, Y , and Z are position (in meters) of the satellites relative
to the center of Earth in rectilinear coordinates. U , V , and W are the velocity components
(in meters per second) in the x, y and z directions, respectively. Earth is assumed to be
stationary in this project.
The parameters in the equations 1 above are given as follows:
• Radius of Earth: Re = 6.37× 106 (m)
• Mass of the Earth: Me = 5.97× 1024 (kg)
• Gravitational constant : G = 6.67408× 10−11 (m3 kg−1 s−2)
• Mass of the satellites: m = 250 (kg)
• Projected area of the satellites: As = 0.25 (m2)
• Air density: ρa = 5.5× 10−12 (kg m−3)
• Drag coefficient: Cd = 2.2
Furthermore, the following quantities will be of interest while analyzing the satellite orbits:
• Altitude of the satellites (m): h =
√
X2 + Y 2 + Z2 −Re
• Speed of the satellites (ms−1): Vmag =
√
U2 + V 2 +W 2
2. Approach
Orbits of numerous satellites will be simulated. The difference among the trajectories
is due to the initial position and velocity of the satellites. The initial position and velocity
components are stored in the text file: satellite data.txt. Download satellite data.txt
from CANVAS.
Using Euler-Cromer method, the equations 1 can be transformed into the following alge-
braic form:
Un+1 = Un −
[
GMe
Xn
(X2n + Y
2
n + Z
2
n)
3/2
+
CdρaAs
2m
Un
√
U2n + V
2
n +W
2
n
]
∆t,
Vn+1 = Vn −
[
GMe
Yn
(X2n + Y
2
n + Z
2
n)
3/2
+
CdρaAs
2m
Vn
√
U2n + V
2
n +W
2
n
]
∆t,
Wn+1 = Wn −
[
GMe
Zn
(X2n + Y
2
n + Z
2
n)
3/2
+
CdρaAs
2m
Wn
√
U2n + V
2
n +W
2
n
]
∆t, (2)
Xn+1 = Xn + Un+1∆t,
Yn+1 = Yn + Vn+1∆t,
Zn+1 = Zn +Wn+1∆t,
where subscript n denotes variables at current time, subscript n+1 denotes variables at time
that is ∆t ahead.
3. Simulation and analysis of satellite trajectories
Perform the following three tasks and put the analysis results inside the project.m script.
Task 1:
Create function satellite.m. This function solves the equations 2 for a given set of
initial conditions. Use ∆t = 1 s. The function should have the following header: function
[T, X, Y, Z, U, V, W] = satellite(Xo, Yo, Zo, Uo, Vo, Wo, Tf) where the inputs
are components of the initial position (Xo, Yo, Zo) and initial velocity (Uo, Vo, Wo), and
the total simulation time (Tf). The outputs are quantities discussed above. All inputs are
scalars while the outputs are vectors.
Use function satellite to obtain the trajectory of the first satellite (ID: 0001) in the
satellite data.txt file for a duration of Tf = 5,000 s. Create figure 1 to plot the trajectory.
Also include the final position of the satellite and Earth surface in each panel. The Earth
surface data is stored in earth topo.mat. A sample script to include the Earth surface is
posted in CANVAS, see the plot earth.m script.
Set: p1a =evalc('help satellite'); p1b ='See figure 1';
Task 2:
Create function read input.m. This function reads the parameters stored in the file
satellite data.txt into MATLAB. The function should have the following declaration:
function [Xo, Yo, Zo, Uo, Vo, Wo] = read input(input filename, sat id) where
input filename is a string denoting the name of the file to be read and sat id is an inte-
ger indicating the satellite ID number. The outputs are the initial position (Xo, Yo, Zo)
and components of initial velocity (Uo, Vo, Wo). When the input sat id is not available
in the file, the function should set all outputs to NaN and display an error warning to screen.
Use function read input and satellite to obtain the trajectories for the first six satel-
lites (ID: 0001-0006) in the satellite data.txt file for a duration of Tf = 12,400 s. Create
figure 2 which includes 2 panels. On the top panel, plot altitude versus time for the six
satellites. Plot speed versus time in the bottom panel. Use function subplot to create the
figure panels. Use different colors to denote the satellites.
Create a 6-element data structure named stat with the following fields:
• sat id: to include a number ranging from 1 to 6 to indicate the satellite ID.
• final position: to include a 3-element vector indicating the final position of the satel-
lites.
• final velocity: to include a 3-element vector indicating the velocity components (U,
V, W) at final time.
• time lmax altitude: to include the times at which the satellites are farthest from
Earth (i.e., times at which the altitude h shows a local maximum).
• orbital period: to include the time that the satellites take to revolve around Earth
for one revolution. Hint: Subtract the times of the first two local maxima of altitude h.
Create a file named report.txt to report the following:
• Your name on the first line
• Your PID on the second line
• A string satellite ID, travel distance (m), orbital period (s) on the third line
• Corresponding values of satellite ID, travel distance and orbital period for each of the
six satellites from line fourth to line nineth. Use single digit for the satellite ID and
format %15.9e for others. Use comma to separate values on the same row.
Set:
p2a =evalc('help read input'); p2b ='See figure 2';
p2c = stat(1); p2d = stat(2); p2e = stat(3);
p2f = stat(4); p2g = stat(5); p2h = stat(6);
p2i = evalc('type report.txt');
Task 3:
Use function read input and satellite to obtain the trajectories for all satellites (ID:
0001-0520) in the satellite data.txt file for a duration of Tf = 6,000 s. Create figure 3 which
includes 2 side-by-side panels:
• Plot the initial position of all the satellites along with the earth surface on the left
panel. Use magenta filled circles to mark the position of the satellites. Use blue filled
circles to mark the position of the satellites that are within 2,000 km distance from
the location at (x = -5.5 ×106 m, y = -3.9 ×106 m, z = 0 m).
• Plot the final position of all the satellites along with the earth surface on the right
panel. Use magenta filled circles to mark the position of the satellites. Use blue filled
circles to mark the final location of satellites that are initially marked in blue in the
left panel.
Set: p3a ='See figure 3'; and answer the following questions using figure 2:
Q1. The satellites move fastest when they are closest or farthest from Earth? Put the
answer in string p3b = '...'.
Q2. As the satellites travel away from earth, does their speed increase or decrease? Put
the answer in string p3c = '...'.
4. Submission instructions
Follow the homework solution template. Remember to clear all, close all, clc,
and fill in your name and PID. Set hw num = 'project'. Create a zip archive named
project.zip. The zip archive should include the following files: project.m, satellite.m,
read input.m, satellite data.txt, report.txt, earth topo.mat and other scripts / func-
tions that you have written for the project. Make sure that you include all necessary
files so that your project.m will run properly. Submit project.zip in CANVAS before
10 PM on Sunday 12/4/2022. Use double precision unless otherwise stated.
Be sure to give all functions a description. All figures need to include axis labels with
correct unit, title and legend (as necessary). Three-dimensional figures need to be set in 3-D
view. The total execution time of your project must be less than 120 s. Zero credit will be
given when your code exceeds this limit.