      function [geog4] = geog4(arg)
% input:
%   arg    = geocentric colatitude (radians)
% output:
%   geog4 = geographic colatitude (radians)
% (n.b. fac=(1-f)**2)
%
      pi2 = pi / 2.0;
      fac = 0.993305621334896;
      geog4 = pi2-atan(cos(arg)/(fac*max(1.e-30,sin(arg))));