function V = orbitalvelocity(R, G, H, units)

% CC Violeta Ivanova, MIT Information Services and Technology
% ORBITALVELOCITY computes the orbital velocity V of a rocket that moves 
% on a circular orbit at altitude H above a planet with mean radius Re and
% gravitational constant g0. 
% V = ORBITALVELOCITY(R, G, H) returns the velocity V computed with the
% formula V = sqrt (G*R^2 / (R+H)), assuming metric units as default.
% V = ORBITALVELOCITY (R, G, H, UNITS) computes the velocity V in English 
% units if the argument UNITS is 'e' or 'English' (not case-sensitive); in 
% metric units if UNITS is 'm' or 'metric'; otherwise it produces an error.


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

if nargin < 3 || nargin > 4
    error('Wrong number of input arguments');
    
elseif nargin == 3

    disp('The solution assumes metric units: planet radius in km, gravity in m/sec^2.')
    V = sqrt(G * R^2 / (R + H) * ucm);
    strV = num2str(V); strH = num2str(H);
    str = ['The orbital velocity for altitude ' strH ' km is: ' strV 'km/hr.'];
    disp(str)
    
else 
    units = lower(units);
     
    switch units
        case {'m', 'metric'}
            V = sqrt(G * R^2 / (R + H) * ucm);
            disp('In metric units: planet radius in km, gravity in m/sec^2.')
            strV = num2str(V); strH = num2str(H);
            str = ['The orbital velocity for altitude ' strH ' km is ' strV ' km/hr.'];
            disp(str)
    
        case {'e', 'english'}
            V = sqrt(G * R^2 / (R + H) * uce);
            strV = num2str(V); strH = num2str(H);
            disp('In English units: planet radius in miles, gravity in ft/sec^2.')
            str = ['The orbital velocity for altitude ' strH ' miles is ' strV ' mi/hr.'];
            disp(str)
        
        otherwise
            error('Unsupported units.')
    end
    
end