MA/CS 321
Fall 2019
1 Goal
Programming Project 1
Due: 11:00am, Monday, Oct 7th. NO late submissions are allowed.
In orbital mechanics, the position of an orbiting body, e.g., a satellite, is determined by solving the Kepler’s equation
M (t) = E (t) − e sin E (t), (1) which relates the eccentric anomaly E, with the eccentricity of the orbit e and the average anomaly M. For
validation purposes, astronomers want to build a library with the routines to find roots of the Kepler’s equation.
2 Detailed description 2.1 Mathematical Model
The position of the satellite is given by the true anomaly ν:
x(t) = r cos ν(t)
where
The true anomaly ν relates to (1) through the equation
e−cosE
(2) (3)
y(t) = r sin ν(t) a(1 − e2) 3 T 2
5 km3 r= 1+ecosν, a= μ 2π , μ=3.986012×10 s2 .
ν=arccos ecosE−1 . (4) and the dependence on time is established through M by
M = 2π t , (5) T
where T is the period of the orbit. To determine the ν, take into account that ν ∈ [0,π] if E ∈ [0,π], and ν ∈ [π,2π] if E ∈ [π,2π].
Then, for all i ∈ {0, . . . , n} find M (ti ) from (5), you have to solve the Kepler equation (1), and finally find
ν(ti) using (4). To solve the Kepler equation, you shall combine Newton’s method and bisection method; a
possible way of doing that is described in Algorithm 1. To sum up, we find the values of all variables via the
order: ti → M(ti) → E(ti) → ν(ti) −→ (x(ti),y(ti)), by plotting (x(ti),y(ti))’s we can get the orbits.
2.2 Routines
You have to code: 1. A routine
function [root, a, b] = bis(a,b,tolb,func,pars)
to find a root of a function f(x) ∈ C(R) using bisection method. We use pars to store the constants M(t)
and e in (1), which can be a vector or other types of array. 2. A routine
(5) (1) (4) (2)(3)
1
MA/CS 321 Fall 2019
Algorithm 1: Root finding algorithm combining bisection method and Newton’s method
INPUT : endpoints of the interval containing the zero a, b; tolerances for the bisection method and
Newton’s method tolb, toln; maximum number of iterations for the Newton’s method
nmax
OUTPUT: eccentric anomaly E
Step 1: Set error = |b − a|; Step 2: while error > toln do
Step 3: set {a, b} = bisection(a, b, tolb);
Step 4: set {E, error} = newton(a + (b − a)/2, toln, nmax); Step 5: if error > toln then set tolb = tolb/2;
else OUTPUT(E); STOP.;
end
OUTPUT(‘Method failed to achieve error ≤ toln’); STOP.
function [root, found, iter] = nwt(seed,toln,nmax,func,dfunc,pars)
to find a root of a function f(x) ∈ C2(R) using Newton’s method. Here found is a logical variable indicating
whether the root is found, and iter is the number of iterations actually used. 3. A routine
function [root, iter, itern] = bisnwt(a,b,tolb,toln,nmax,func,dfunc,pars)
to find a root of a function f(x) ∈ C2(R) using Algorithm 1. Here iter is the number of iterations used, and
itern is the number of iterations used in the Newton’s method. 4. An application
function [orbit] = tanom(T,e,n,varargin)
that solves the Kepler equation and gives the position of the satellite in its orbital plane. This function expects the arguments:
• T: the period of the orbit in hours.
• e: the eccentricity of the orbit.
• n: the number of points where we want to localize the satellite on its orbit.
Moreover, the following arguments can optionally be passed to the function through varargin. (see math- works):
• tolb: the error for bisection method;
• toln: the error for Newton method;
• nmax: the maximum iterations that we allow to Newton method.
The output argument of the program, orbit, must be a matrix of the format
t0 E0 ν0 r0 x0 y0 t1 E1 ν1 r1 x1 y1 . . . . . .
tn En νn rn xn yn
where ti = iT/n, Ei is the solution to Kepler Equation, νi is the true anomaly corresponding to ti, ri is the
distance between the satellite and the earth at ti and (xi,yi) are the cartesian coordinates. Each function must be in a different .m file, with the same name as the function.
2
MA/CS 321 Fall 2019
3 Summary
You are required to submit to “Project 1” in Canvas a .zip file containing the following files:
1. The source files “bis.m”, “nwt.m”, “bisnwt.m”, “tanom.m” described above1.
2. A .pdf file with a brief report of the project. The report should include the following sections:
• User Guide: A concise description of all the routines and applications: what do they do, calling syntax, output format, and return values or messages.
• Implementation:
– The argument toln in the routine tanom.m represents the error that we allow in the approximation of the Kepler equation. If we want the error in the coordinates of the orbit to be smaller than 10−6, which value should we use for toln? Justify your answer.
– Intuitively, the root-finding routines are meant for x ∈ R and thus, the points in the orbit can be found one by one, applying Algorithm 1 and the transformations described in Section 2.1 for each time of the period. A slight change, though, would allow us to work with a vector x ∈ Rn and find all the points of the orbit at once. Discuss about the differences between the two approaches.
– Describe the tests that you have done to validate your code, what you have discovered, the problems encountered, and their respective solutions. You must include at least the following test
T =4, e=0.25, n=100, and tolerance 1e-12 in “bisnwt.m”.
• Solutions: Plot the orbits that you have obtained with your selected initial data and parameters.
1If you use python or other programming languages, please create the files/functions with the same names. 3