function [r,g] = georad(m,g);

% GEORAD  Radius of spherical approximation of an ellipsoid
%
% Radius = GEORAD( Mode , [Geoid] )
%
% Mode:  'a'  Authalic (Equal Surface Area)  (default)
%        'v'  Equal Volume
%        'r'  Rectifying (Equal Meridional Distances )
%
% Geoid = [ MajorAxis Eccentricity ],  optional, default: WGS84
%
% [ Radius , Geoid ] = GEORAD( Mode , [Geoid] )
%
%------------------------------------------------------------------
%
% See also: GEOLAT
%


Nin = nargin;

if Nin < 1
   m = '';
end

if Nin < 2
   g = [];
end

if isempty(m)
   m = 'a';
end

%******************************************************************
if ~isempty(g)
%******************************************************************
    if ~( ( ndims(g) == 2 ) & ( size(g,2) == 2 ) );
        error('Geoid must have 2 Columns: [ MajorAxis Eccentricity ].')
    end
    a = g(:,1);
    e = g(:,2);
%******************************************************************
else
%******************************************************************
% WGS84

a = 6378137.;
b = 6356752.314;        %   A * sqrt(1-E*E);
e = sqrt(a*a - b*b)/a;  %   E = 0.081819191;

g = [ a  e ];

%******************************************************************
end
%******************************************************************

if e == 0
   r = a;
   return
end

switch upper(m(1))

  %****************************************************************
  case 'A'  % Authalic
  %****************************************************************

    r = a/2 .* sqrt( 2 + ( 1./e - e ) .* log((1+e)./(1-e)) );

  %****************************************************************
  case 'V'   % EqualVolume
  %****************************************************************

    f = 1 - sqrt( 1 - e.^2 );  % Flattening

    r = a .* ( 1 - f/3 - f.^2 / 9 );

  %****************************************************************
  case 'R'   % Rectifying
  %****************************************************************

    n = ( 1 - sqrt(1-e.^2) ) ./ ( 1 + sqrt(1-e.^2) ); % (a-b) / (a+b);

    f = ( 1 - n ) .* ( 1 - n.^2 );                    % 8*a*b^2 / (a+b)^3;

    r = a .* f .* ( 1 + 9/4 * n.^2 + 225/64 * n.^4 );

  %****************************************************************
  otherwise
  %****************************************************************

    error('Invalid Mode.');

end

%******************************************************************
