MatLab

Appendix 2 – Week 2 Assisting MatLab Scripts
2A – Kepler’s Universal solution for Universal Anomaly
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function x = kepler_U(dt, ro, vro, a)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%{
This function uses Newton's method to solve the universal
Kepler equation for the universal anomaly.
mu - gravitational parameter (km^3/s^2)
x - the universal anomaly (km^0.5)
dt - time since x = 0 (s)
ro - radial position (km) when x = 0
vro - radial velocity (km/s) when x = 0
a - reciprocal of the semimajor axis (1/km)
z - auxiliary variable (z = a*x^2)
C - value of Stumpff function C(z)
S - value of Stumpff function S(z)
n - number of iterations for convergence
nMax - maximum allowable number of iterations
User M-functions required: stumpC, stumpS
% ----------------------------------------------
global mu
%...Set an error tolerance and a limit on the number of iterations:
error = 1.e-8;
nMax = 1000;
%...Starting value for x:
x = sqrt(mu)*abs(a)*dt;
%...Iterate on Equation 3.65 until until convergence occurs within
%...the error tolerance:
n = 0;
ratio = 1;
while abs(ratio) > error && n <= nMax
    n = n + 1;
    C = stumpC(a*x^2);
    S = stumpS(a*x^2);
    F = ro*vro/sqrt(mu)*x^2*C + (1 - a*ro)*x^3*S + ro*x - sqrt(mu)*dt;
    dFdx = ro*vro/sqrt(mu)*x*(1 - a*x^2*S) + (1 - a*ro)*x^2*C + ro;
    ratio = F/dFdx;
    x = x - ratio;
end
%...Deliver a value for x, but report that nMax was reached:
if n > nMax
    fprintf('\n **No. iterations of Kepler''s equation = %g', n)
    fprintf('\n F/dFdx = %g\n', F/dFdx)
end
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2B – Computation of f and g Lagrange Coefficients
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function [f, g] = f_and_g(x, t, ro, a)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%{
This function calculates the Lagrange f and g coefficients.
mu - the gravitational parameter (km^3/s^2)
a - reciprocal of the semimajor axis (1/km)
ro - the radial position at time to (km)
t - the time elapsed since ro (s)
x - the universal anomaly after time t (km^0.5)
f - the Lagrange f coefficient (dimensionless)
g - the Lagrange g coefficient (s)
User M-functions required: stumpC, stumpS
%}
% ----------------------------------------------
global mu
z = a*x^2;
%...Equation 3.69a:
f = 1 - x^2/ro*stumpC(z);
%...Equation 3.69b:
g = t - 1/sqrt(mu)*x^3*stumpS(z);
end
2C – Computation of the derivation of f and g Lagrange Coefficients
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function [fdot, gdot] = fDot_and_gDot(x, r, ro, a)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%{
This function calculates the time derivatives of the
Lagrange f and g coefficients.
mu - the gravitational parameter (km^3/s^2)
a - reciprocal of the semimajor axis (1/km)
ro - the radial position at time to (km)
t - the time elapsed since initial state vector (s)
r - the radial position after time t (km)
x - the universal anomaly after time t (km^0.5)
fdot - time derivative of the Lagrange f coefficient (1/s)
gdot - time derivative of the Lagrange g coefficient (dimensionless)
User M-functions required: stumpC, stumpS
%}
% --------------------------------------------------
global mu
z = a*x^2;
%...Equation 3.69c:
fdot = sqrt(mu)/r/ro*(z*stumpS(z) - 1)*x;
%...Equation 3.69d:
gdot = 1 - x^2/r*stumpC(z);
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2D – Computation of the Stumpff Function C(z)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function c = stumpC(z)
% ~~~~~~~~~~~~~~~~~~~~~~
%{
This function evaluates the Stumpff function C(z) according
to Equation 3.53.
z - input argument
c - value of C(z)
User M-functions required: none
%}
% ----------------------------------------------
if z > 0
    c = (1 - cos(sqrt(z)))/z;
elseif z < 0
    c = (cosh(sqrt(-z)) - 1)/(-z);
else
    c = 1/2;
end
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2E – Computation of the Stumpff Function S(z)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
function s = stumpS(z)
% ~~~~~~~~~~~~~~~~~~~~~~
%{
This function evaluates the Stumpff function S(z) according
to Equation 3.52.
z - input argument
s - value of S(z)
User M-functions required: none
%}
% ----------------------------------------------
if z > 0
    s = (sqrt(z) - sin(sqrt(z)))/(sqrt(z))^3;
elseif z < 0
    s = (sinh(sqrt(-z)) - sqrt(-z))/(sqrt(-z))^3;
else
    s = 1/6;
end
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~