function [NORTH, EAST] = dms2utm(lad, lam, las, lod, lom, los)
%[NORTH, EAST] = dms2utm(lad, lam, las, lod, lom, los)
%	lad, lam, las  are latitude in degrees, minutes, seconds north
%	   (set lad negative for southern hemisphere)
%	lod, lom, los  are longitude in degrees, minutes, seconds west
%	   (set lod negative or > 179 for eastern hemisphere)
%	These parameters may be arrays of equal dimension.
%       Output is in NAD27 meters north and east of zone reference lines.
%       Formulas adapted from "UTMS" by E. Carlson and T. Vincenty
%	   (ftp://ftp.ngs.noaa.gov/pub/pcsoft/utms/ )

%% WGS84 / NAD83 Ellipsoid data
% er = 6378137;
% rf = 298.257222101;

%% NAD27 Ellipsoid data
er = 6378206.4;
rf = 294.978698;

rad = 180/pi;

fe = 500000;   % false easting moves reference line west of zone meridian

sf = 0.9996;   % scale factor for UTM conversion

f = 1 / rf;
esq = f + f - f^2;
eps = esq / (1 - esq);

phi = (abs(lad) + (lam + las / 60) / 60) / rad;
fn = (lad < 0) * 10000000;	% false northing in southern hemisphere

ts = tan(phi).^2;
ets = eps * cos(phi).^2;

theta = (lod + (lom + los / 60) / 60) / rad;

western = (mod(lod,360) < 180);	% western hemisphere
zone = western .* floor(30 - mod(lod,360) / 6) + ...
       (1 - western) .* floor(90 - mod(lod,360) / 6);
mid = western .* (183 - 6 * zone) / rad + ...
      (1 - western) .* (543 - 6 * zone) / rad;   % longitude of zone meridian

l = (theta - mid) .* cos(phi);
ls = l.^2;

pr = (1 - f) * er;
en = (er - pr) / (er + pr);

r = er * (1 - en) * (1 - en^2) * (1 + 2.25 * en^2 + (225 / 64) * en^4);

a = -1.5 * en + (9 / 16) * en^3;
b = 0.9375 * en^2 - (15 / 32) * en^4;
c = - (35 / 48) * en^3;

omega = phi + a * sin(2 * phi) + b * sin(4 * phi) + c * sin(6 * phi);
s = r * omega * sf;

rn = sf * er ./ sqrt(1 - esq * sin(phi).^2);

a1 = -rn;
a2 = rn .* tan(phi) / 2;
a3 = (1 - ts + ets) / 6;
a4 = (5 - ts + ets .* (9 + 4 * ets)) / 12;
a5 = (5 + ts .* (ts - 18) + ets .* (14 - 58 * ts)) / 120;
a6 = (61 + ts .* (ts - 58) + ets .* (270 + 330 * ts)) / 360;
a7 = (61 - 479 * ts + 179 * ts.^2 - ts.^3) / 5040;

NORTH = abs(s + a2 .* ls .* (1 + ls .* (a4 + a6 .* ls)) - fn);
EAST = fe + a1 .* l .* (1 + ls .* (a3 + ls .* (a5 + a7 .* ls)));