function  [x,y,utmzone] = wgs2utm(Lat,Lon)
% This function converts the vectors of Lat/Lon coordinates to those of UTM.
% Inputs:
%    Lat (WGS84 Latitude vector)  in decimal degrees 
%    Lon (WGS84 Longitude vector) in decimal degrees
% Outputs:
%    x       - UTM easting in meters
%    y       - UTM northing in meters
%    utmzone - UTM longitudinal zone
%
% Example:
%    Lat=[48.866667; 34.05;   -36.85];
%    Lon=[2.333056;  -118.25; 174.783333];
%    [x,y,utmzone]=wgs2utm(Lat,Lon)
%       returns
% x =
%   1.0e+005 *
%     4.5109
%     3.8463
%     3.0237
% y =
%   1.0e+006 *
%     5.4128
%     3.7684
%     5.9195
% utmzone =
% 31U
% 11S
% 60H
%
% Source: DMA Technical Manual 8358.2, Fairfax, VA

%% Argument checking
error(nargchk(2,2,nargin));         % 2 arguments are required
n1=size(Lat);
n2=size(Lon);
if (n1~=n2)
   error('Lat and Lon should have same size'); return
end

%% Converting coordinates to radians
lat = Lat*pi/180;
lon = Lon*pi/180;

%% WGS84 parameters
a = 6378137;            % semi-major axis of the Earth ellipsoid
b = 6356752.314245;     % semi-minor axis of the Earth ellipsoid
e = sqrt(1-(b/a)^2);    % first eccentricity

%% UTM parameters
Lon0 = floor(Lon/6)*6+3;    % reference longitude in degrees
lon0 = Lon0.*pi/180;        % reference longitude in radians
k0 = 0.9996;                % scale on central meridian

FE = 500000;                % false easting
FN = (Lat < 0).*10000000;   % false northing 

%% Equations parameters
eps = e^2/(1-e^2);          % squared second eccentricity
% N is the radius of curvature of the earth perpendicular to meridian plane
% Also, distance from a point to polar axis
N = a./sqrt(1-e^2*sin(lat).^2);
T = tan(lat).^2;             
C = eps*(cos(lat)).^2;
A = (lon-lon0).*cos(lat);                       
% M: true distance along the central meridian from the equator to lat
M = a*(  (1 - e^2/4 - 3*e^4/64 - 5*e^6/256)* lat         ...
        -(3*e^2/8 + 3*e^4/32 + 45*e^6/1024)* sin(2*lat) ...
        +(15*e^4/256 + 45*e^6/1024)        * sin(4*lat) ...
        -(35*e^6/3072 )                    * sin(6*lat));

%% Computing easting
x = FE + k0*N.*(                             A      ...
               +(1-T+C)                   .* A.^3/6 ...
               +(5-18*T+T.^2+72*C-58*eps) .* A.^5/120);
               
%% Computing northing 
y = FN + k0*M + k0*N.*tan(lat).*(                                 A.^2/2  ...
                                 +(5-T+9*C+4*C.^2)             .* A.^4/24 ...
                                 +(61-58*T+T.^2+600*C-330*eps) .* A.^6/720);
                                 
%% Determining the UTM zone
lonzone = floor(Lon0./6)+31;
latband = floor((Lat+80)./8)+1;
LB='CDEFGHJKLMNPQRSTUVWX';
latband=LB(latband)';
utmzone = [num2str(lonzone,'%02g') latband];