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

% ORBITALVELOCITY computes a rocket's circular orbital velocity 
% at altitude H above a planet with mean radius R and gravity G. 
% V = ORBITALVELOCITY(R, G, H) returns the velocity V computed with J.
% Kepler's formula V = sqrt (G*R^2/(R+H)), assuming metric units.
% 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.
% CC Violeta Ivanova, violeta@mit.edu

% 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/h.'];
    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/h.'];
            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 ' mph.'];
            disp(str)
        
        otherwise
            error('Unsupported units.')
    end
    
end