function geo = utm2geo( utm, cm, fe, phi_0, fn, a, finv )

%UTM2GEO      Converts UTM coordinates to geodetic coordinates.
%   GEO = UTM2GEO( UTM, CM, FE, PHI_0, FN, A, FINV ) converts M by 2
%   matrix of UTM coordinates [northing easting; ...] in meters into
%   geodetic coordinates [latitude longitude; ...] in degrees, North
%   and East positive.  CM and FE are the central meridian of the
%   desired UTM Zone and the false easting of the projection (default
%   Zone 19 at W69 deg., i.e. CM = -69, FE = 500,000 m).  PHI_0 and FN
%   are the latitude of origin and false northing of the projection
%   (default PHI_0 = 0, FN = 0 m).  If unspecified, the semi-major
%   axis and inverse flattening default to the WGS84 survey parameters
%   (A = 6378137 and FINV = 298.257223563).  This function calls ELLIPARC,
%   FOOTPOINT, and PRIMERAD.

%----------------------------------------------------------------------
%			  Copyright (C) 1995
%		Massachussetts Institute of Technology
%	       Sea Grant Underwater Vehicles Laboratory
%----------------------------------------------------------------------
%     $Revision: 1.1.1.1 $			$State: Exp $
%     $Author: cvsadmin $		$Date: 2007-07-13 03:02:12 $
%----------------------------------------------------------------------

if nargin < 7, finv = 298.257223563; end
if nargin < 6, a = 6378137; end
if nargin < 5, fn = 0; end
if nargin < 4, phi_0 = 0; end
if nargin < 3, fe = 500000; end
if nargin < 2, cm = -69; end

global Y k0;
k0 = 0.9996;			% scale factor at central meridian

phi_f = zeros(size(utm,1),1);
				% solve for footpoint latitude phi_f
for i = 1:length(phi_f)
    Y = utm(i,1) - fn + k0*elliparc( phi_0*pi/180, a, finv );
    phi_f(i) = fzero( 'footpoint', phi_f(i), 1e-10 );
end

S = elliparc( phi_f, a, finv );	% length along elliptical arc
N = primerad( phi_f, a, finv );	% radius prime vertical section

t = tan( phi_f );		% intermediate variables
t2 = t.^2;
cp_f = cos( phi_f );
h2 = (2*finv -1)/(finv - 1)^2*cp_f.^2;
x_n = (utm(:,2) - fe)./(k0*N);
x_n2 = x_n.^2;

geo = zeros(size(utm));         % preallocate geo

                                % latitude
geo(:,1) = phi_f ...
         - t.*x_n2.*((1 + h2)/2 ...
                   + x_n2.*((5 + 3*t2 + 6*h2.*(1 - t2) ...
                           - 3*h2.*h2.*(1 + 3*t2))/24 ...
                           - x_n2/720.*(61 + t2.*(90 + 45*t2) + 107*h2 ...
                                      - 162*t2.*h2 - 45*t2.*t2.*h2)));
                                % longitude
geo(:,2) = x_n.*(1 - x_n2.*((1 + 2*t2 + h2)/6 ...
               - x_n2.*(5 + t2.*(28 + 24*t2) + h2.*(6 + 8*t2))/120))./cp_f;


geo = geo*180/pi;		% convert to degrees
geo(:,2) = geo(:,2) + cm;	% add central meridian
