function  [x,y,utmzone] = deg2utm(Lat,Lon)

 %Modified by : Marc Bosch (speed up the code by vectorizing operations) JHUAPL

% -------------------------------------------------------------------------
% [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Description: Function to convert lat/lon vectors into UTM coordinates (WGS84).
% Some code has been extracted from UTM.m function by Gabriel Ruiz Martinez.
%
% Inputs:
%    Lat: Latitude vector.   Degrees.  +ddd.ddddd  WGS84
%    Lon: Longitude vector.  Degrees.  +ddd.ddddd  WGS84
%
% Outputs:
%    x, y , utmzone.   See example
%
% Example 1:
%    Lat=[40.3154333; 46.283900; 37.577833; 28.645650; 38.855550; 25.061783];
%    Lon=[-3.4857166; 7.8012333; -119.95525; -17.759533; -94.7990166; 121.640266];
%    [x,y,utmzone] = deg2utm(Lat,Lon);
%    fprintf('%7.0f ',x)
%       458731  407653  239027  230253  343898  362850
%    fprintf('%7.0f ',y)
%      4462881 5126290 4163083 3171843 4302285 2772478
%    utmzone =
%       30 T
%       32 T
%       11 S
%       28 R
%       15 S
%       51 R
%
% Example 2: If you have Lat/Lon coordinates in Degrees, Minutes and Seconds
%    LatDMS=[40 18 55.56; 46 17 2.04];
%    LonDMS=[-3 29  8.58;  7 48 4.44];
%    Lat=dms2deg(mat2dms(LatDMS)); %convert into degrees
%    Lon=dms2deg(mat2dms(LonDMS)); %convert into degrees
%    [x,y,utmzone] = deg2utm(Lat,Lon)
%
% Author: 
%   Rafael Palacios
%   Universidad Pontificia Comillas
%   Madrid, Spain
% Version: Apr/06, Jun/06, Aug/06, Aug/06
% Aug/06: fixed a problem (found by Rodolphe Dewarrat) related to southern 
%    hemisphere coordinates. 
% Aug/06: corrected m-Lint warnings
%-------------------------------------------------------------------------


error(nargchk(2, 2, nargin));  %2 arguments required
n1=length(Lat);
n2=length(Lon);
if (n1~=n2)
   error('Lat and Lon vectors should have the same length');
end


% Memory pre-allocation
%
x=zeros(n1,1);
y=zeros(n1,1);
utmzone(n1,:)='60 X';

 la=Lat;
    if (la<-72), Letra='C';
   elseif (la<-64), Letra='D';
   elseif (la<-56), Letra='E';
   elseif (la<-48), Letra='F';
   elseif (la<-40), Letra='G';
   elseif (la<-32), Letra='H';
   elseif (la<-24), Letra='J';
   elseif (la<-16), Letra='K';
   elseif (la<-8), Letra='L';
   elseif (la<0), Letra='M';
   elseif (la<8), Letra='N';
   elseif (la<16), Letra='P';
   elseif (la<24), Letra='Q';
   elseif (la<32), Letra='R';
   elseif (la<40), Letra='S';
   elseif (la<48), Letra='T';
   elseif (la<56), Letra='U';
   elseif (la<64), Letra='V';
   elseif (la<72), Letra='W';
   else Letra='X';
   end
   lo=Lon;

   sa = 6378137.000000 ; sb = 6356752.314245;
         
   %e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
   e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
   e2cuadrada = e2 ^ 2;
   c = ( sa ^ 2 ) / sb;
   %alpha = ( sa - sb ) / sa;             %f
   %ablandamiento = 1 / alpha;   % 1/f

   lat = la * ( pi / 180 );
   lon = lo * ( pi / 180 );

   Huso = fix( ( lo / 6 ) + 31);
   S = ( ( Huso * 6 ) - 183 );
   deltaS = lon - ( S * ( pi / 180 ) );
    
   a = cos(lat) .* sin(deltaS);
   epsilon = 0.5 * log( ( 1 +  a) ./ ( 1 - a ) );
   nu = atan( tan(lat) ./ cos(deltaS) ) - lat;
   v = ( c ./ ( ( 1 + ( e2cuadrada .* ( cos(lat) ) .^ 2 ) ) ) .^ 0.5 ) * 0.9996;
   ta = ( e2cuadrada / 2 ) .* epsilon .^ 2 .* ( cos(lat) ) .^ 2;
   a1 = sin( 2 * lat );
   a2 = a1 .* ( cos(lat) ) .^ 2;
   j2 = lat + ( a1 / 2 );
   j4 = ( ( 3 * j2 ) + a2 ) / 4;
   j6 = ( ( 5 * j4 ) + ( a2 .* ( cos(lat) ) .^ 2) ) / 3;
   alfa = ( 3 / 4 ) * e2cuadrada;
   beta = ( 5 / 3 ) * alfa .^ 2;
   gama = ( 35 / 27 ) * alfa .^ 3;
   Bm = 0.9996 * c .* ( lat - alfa .* j2 + beta .* j4 - gama .* j6 );
   xx = epsilon .* v .* ( 1 + ( ta / 3 ) ) + 500000;
   yy = nu .* v .* ( 1 + ta ) + Bm;
   pos=find(yy<0);
   yy(pos) = 9999999+yy(pos);
   if (yy<0)
       yy=9999999+yy;
   end
   x=xx;y=yy;
   utmzone=sprintf('%02d %c',Huso,Letra);
% 
% 
% % Main Loop
% %
% for i=1:n1
%    la=Lat(i);
%    lo=Lon(i);
% 
%    sa = 6378137.000000 ; sb = 6356752.314245;
%          
%    %e = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sa;
%    e2 = ( ( ( sa ^ 2 ) - ( sb ^ 2 ) ) ^ 0.5 ) / sb;
%    e2cuadrada = e2 ^ 2;
%    c = ( sa ^ 2 ) / sb;
%    %alpha = ( sa - sb ) / sa;             %f
%    %ablandamiento = 1 / alpha;   % 1/f
% 
%    lat = la * ( pi / 180 );
%    lon = lo * ( pi / 180 );
% 
%    Huso = fix( ( lo / 6 ) + 31);
%    S = ( ( Huso * 6 ) - 183 );
%    deltaS = lon - ( S * ( pi / 180 ) );
% 
%    if (la<-72), Letra='C';
%    elseif (la<-64), Letra='D';
%    elseif (la<-56), Letra='E';
%    elseif (la<-48), Letra='F';
%    elseif (la<-40), Letra='G';
%    elseif (la<-32), Letra='H';
%    elseif (la<-24), Letra='J';
%    elseif (la<-16), Letra='K';
%    elseif (la<-8), Letra='L';
%    elseif (la<0), Letra='M';
%    elseif (la<8), Letra='N';
%    elseif (la<16), Letra='P';
%    elseif (la<24), Letra='Q';
%    elseif (la<32), Letra='R';
%    elseif (la<40), Letra='S';
%    elseif (la<48), Letra='T';
%    elseif (la<56), Letra='U';
%    elseif (la<64), Letra='V';
%    elseif (la<72), Letra='W';
%    else Letra='X';
%    end
% 
%    a = cos(lat) * sin(deltaS);
%    epsilon = 0.5 * log( ( 1 +  a) / ( 1 - a ) );
%    nu = atan( tan(lat) / cos(deltaS) ) - lat;
%    v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
%    ta = ( e2cuadrada / 2 ) * epsilon ^ 2 * ( cos(lat) ) ^ 2;
%    a1 = sin( 2 * lat );
%    a2 = a1 * ( cos(lat) ) ^ 2;
%    j2 = lat + ( a1 / 2 );
%    j4 = ( ( 3 * j2 ) + a2 ) / 4;
%    j6 = ( ( 5 * j4 ) + ( a2 * ( cos(lat) ) ^ 2) ) / 3;
%    alfa = ( 3 / 4 ) * e2cuadrada;
%    beta = ( 5 / 3 ) * alfa ^ 2;
%    gama = ( 35 / 27 ) * alfa ^ 3;
%    Bm = 0.9996 * c * ( lat - alfa * j2 + beta * j4 - gama * j6 );
%    xx = epsilon * v * ( 1 + ( ta / 3 ) ) + 500000;
%    yy = nu * v * ( 1 + ta ) + Bm;
% 
%    if (yy<0)
%        yy=9999999+yy;
%    end
% 
%    x(i)=xx;
%    y(i)=yy;
%    utmzone(i,:)=sprintf('%02d %c',Huso,Letra);
% end
