% IAP 2006 Introduction to MATLAB
% Interface and Basics

% EXAMPLE 2 Orbital Velocity
% Based on http://exploration.grc.nasa.gov/education/rocket/rktrflght.html

% A. Matrices and vectors 
% Compute a rocket's velocity V for a circular orbit around a planet.
% Re is the planet's mean radius.
% g0 is the planet's surface gravitational constant. 
% h is the altitude of the circular orbit.
% V can be computed using the formula developed by Johannes Kepler:
% V = sqrt(g0*Re^2/(Re+h))

% B. Compute V for the Earth, Moon, and Mars.

% SOLUTION

% B1. Create a matrix with the values of Re for the Earth, Moon, and Mars.
% The first column is in English units (miles), the second column is in
% metric units (km).
Re = [3963 1079 2111; 6376 1736 3396]'

% B2. Create a matrix with the values of g0 for the Earth, Moon, and Mars.
% The first column is in English units (ft/sec^2).
% The second column is in metric units (m/sec^2).
g0 = [32.2 9.814; 5.3 1.615; 12.1 3.688]

% B3. Create coefficients for unit conversion.
% For English units:
uce = 3600^2/5280;
% For metric units:
ucm = 3600^2/1000;

% B4. Solve for orbital altitude h = 100 miles = 160 km.

% In English units for the Earth, Moon, and Mars:
g0e = g0(:, 1)
Ree = Re(:, 1)
Ve = sqrt(g0e .* Ree .^2 ./ (Ree+100)*uce)

% In metric units for the Earth, Moon, and Mars:
g0m = g0(:, 2)
Rem = Re(:, 2)
Vm = sqrt(g0m .* Rem .^2 ./ (Rem+160)*ucm)

% B5. Combine solutions for V in a matrix with rows for the Earth, Moon, and Mars.
% The first column is in English units (miles/hr).
% The second column is in metric units (km/hr).
V = [Ve Vm]

% C. Polynomials 
% Solve a polynomial equation to find V for altitude h=250 miles above Earth.
% The formula for orbital velocity can be rewritten as:
% V^2 - g0*Re^2/(Re+h) = 0
% This is a polynomial equation for V of the type:
% p1*V^2 + p2*V + p3 = 0

% SOLUTION

% C1. Create a vector with the coefficients of the polynomial equation:
h=250; 
p3 = -g0(1,1)*Re(1,1)^2*uce/(Re(1,1)+h)
p = [1 0 p3]
% where g0(1,1) and Re(1,1) are the first elements in the matrices g0 and Re, respectively. 
% and uce is the unit conversion multiplier for English units.

% C2. Use the built-in function roots to compute the velocity V in mi/hr for h=250 miles
v = roots(p)

% C3. Check the solution using the built-in function polyval
P = polyval(p, v)
% The vector P should have elements that are equal to or very close to zero.