function [out1,out2,savepts] = tranmerc(mstruct,in1,in2,object,direction,savepts)

%TRANMERC  Transverse Mercator Projection
%
%  This conformal projection is the transverse form of the Mercator
%  projection and is also known as the Gauss-Krueger pojection.  It is not
%  equal area, equidistant, or perspective.
%
%  The scale is constant along the central meridian, and increases to the
%  east and west. The scale at the entral meridian can be set true to
%  scale, or reduced slightly to render the mean scale of the overall map
%  more nearly correct.
%
%  The uniformity of scale along its centeral meridian makes Transverse
%  Mercator an excellent choice for mapping areas that are elongated
%  north-to-south.  Its best known application is the definition of
%  Universal Transverse Mercator (UTM) coordinates.  Each UTM zone spans
%  only 6 degrees of longitude, but the northern half extends from the
%  equator all the way to 84 degrees north and the southern half extends
%  from 80 degrees south to the equator. Other map grids based on
%  Transverse Mercator include many of the state plane zones in the U.S.,
%  and the U.K. National Grid.

%   Copyright 2002-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $ $Date: 2004/12/18 07:46:49 $

if nargin == 1                  %  Set the default structure entries
	  mstruct.trimlat = angledim([-80 80],'degrees',mstruct.angleunits);
      mstruct.trimlon = angledim([-20 20],'degrees',mstruct.angleunits);
	  mstruct.mapparallels = [];
      mstruct.nparallels   = 0;
	  mstruct.fixedorient  = [];
	  mstruct.scalefactor  = 1;
      mstruct.aspect = 'normal';
      out1 = mstruct;
      return;
elseif nargin ~= 5 & nargin ~= 6
      error('Incorrect number of arguments')
end

if ~strmatch(mstruct.aspect,'normal')
    warning('Aspect setting ''%s'' ignored by Transverse Mercator projection.',...
            mstruct.aspect);
end

%  Extract the necessary projection data and convert to radians
units     = mstruct.angleunits;
origin    = angledim(mstruct.origin,   units,'radians');
trimlat   = angledim(mstruct.flatlimit,units,'radians');
trimlon   = angledim(mstruct.flonlimit,units,'radians');

switch direction
    case 'forward'
        lat  = angledim(in1,units,'radians');   %  Convert to radians
        long = angledim(in2,units,'radians');
	
        long = long - origin(2);   %  Rotate to new longitude origin
        [lat,long,clipped] = clipdata(lat,long,object);    %  Clip at the date line
        [lat,long,trimmed] = trimdata(lat,trimlat,long,trimlon,object);
	
        %  Construct the structure of altered points
        savepts.trimmed = trimmed;
        savepts.clipped = clipped;
	
        %  Perform the projection calculations
        [x,y] = transmercfwd(lat, long,...
                       mstruct.geoid, mstruct.scalefactor, origin(1), 0,...
                       mstruct.falseeasting, mstruct.falsenorthing);
	
        %  Set the output variables
        out1 = x;
        out2 = y;
	
    case 'inverse'
        x = in1;
        y = in2;

        % Perform inverse projection
        [lat,long] = transmercinv(x, y,...
                            mstruct.geoid, mstruct.scalefactor, origin(1), 0,...
                            mstruct.falseeasting, mstruct.falsenorthing);
	
        %  Undo trims and clips
        [lat,long] = undotrim(lat,long,savepts.trimmed,object);
        [lat,long] = undoclip(lat,long,savepts.clipped,object);
	
        %  Rotate to Greenwhich and transform to desired units
        long = long + origin(2);
        out1 = angledim(lat, 'radians', units);
        out2 = angledim(long,'radians', units);
       
    otherwise
        error('Unrecognized direction string')
end

%  Some operations on NaNs produce NaN + NaNi.  However operations
%  outside the map may product complex results and we don't want
%  to destroy this indicator.

indx = find(isnan(out1) | isnan(out2));
out1(indx) = NaN;
out2(indx) = NaN;
		
%----------------------------------------------------------------------
% Transverse Mercator formulas from the publication:
%   A guide to coordinate systems in Great Britain, Ordnance Survey.
%----------------------------------------------------------------------

function [x,y] = transmercfwd(phi, lambda,...
                              ellipsoid, scalefactor, phi0, lambda0, x0, y0)

F0 = scalefactor;
a = ellipsoid(1);
e2 = (ellipsoid(2))^2;

sinphi = sin(phi);
cosphi = cos(phi);
sinphi2 = sinphi.^2;
cosphi2 = cosphi.^2;

dlam2 = (lambda - lambda0).^ 2;

nu = a * F0 ./ sqrt(1 - e2 * sinphi2);
nuOVERrho = (1 - e2 * sinphi2)/(1 - e2);
eta2 = nuOVERrho - 1;

yc = ycentral(phi, ellipsoid, scalefactor, phi0, y0);

x = x0 + nu .* cosphi .*  (lambda - lambda0) .* (1 ...
    + dlam2 .* (((nuOVERrho) .* cosphi2 - sinphi2)/6 ...
    + dlam2 .* (((5 + 14 * eta2) .* cosphi2.^2 ...
                     - (18 + 58 * eta2) .* cosphi2 .* sinphi2 + sinphi2.^2))/120));

y = yc + nu .* sinphi .* cosphi .* dlam2 .* (1/2 ...
    + dlam2 .* (((5 + 9 * eta2) .* cosphi2 - sinphi2)/24 ...
    + dlam2 .* (61 * cosphi2.^2 ...
                     - 58 * cosphi2 .* sinphi2 + sinphi.^2)/720));

%-----------------------------------------------------------------------

function [phi,lambda] = transmercinv(x, y,...
                                     ellipsoid, scalefactor, phi0, lambda0, x0, y0)

F0 = scalefactor;
a = ellipsoid(1);
e2 = (ellipsoid(2))^2;

phiP = phicentral(y, ellipsoid, scalefactor, phi0, y0);

sinphiP2 = sin(phiP).^2;
nu = a * F0 ./ sqrt(1 - e2 * sinphiP2);
nuOVERrho = (1 - e2 * sinphiP2)/(1 - e2);
eta2 = nuOVERrho - 1;

tanphiP  = tan(phiP);
tanphiP2 = tanphiP.^2;
dx2OVERnu2 = ((x - x0)./nu).^2;

phi = phiP + tanphiP .* nuOVERrho .* dx2OVERnu2 .* (-1/2 ...
      + dx2OVERnu2 .* ((5 + eta2 + tanphiP2 .* (3 - 9 * eta2))/24 ...
      + dx2OVERnu2 .* -(61 + tanphiP2 .* (90 + 45 * tanphiP2))/720));

lambda = lambda0 + sec(phiP) .* ((x - x0)./nu) .* (1 ...
         + dx2OVERnu2 .* (-(nuOVERrho + 2 * tanphiP2)/6 ...
         + dx2OVERnu2 .* ((5 + tanphiP2 .* (28 + 24 * tanphiP2))/120 ...
         + dx2OVERnu2 .* -(61 + tanphiP2 .* (662 + tanphiP2 .* (1320 + 720 * tanphiP2)))/5040)));

%-----------------------------------------------------------------------

function yc = ycentral(phi, ellipsoid, scalefactor, phi0, y0)
% Calculate the northing of the point (phi,lambda0), which lies on the
% central meridian.

F0 = scalefactor;
b = minaxis(ellipsoid);
n = ecc2n(ellipsoid);

dphi = phi - phi0;
pphi = phi + phi0;

yc = y0 + b * F0 ...
          * ( ((4 + n.*(4 + 5*n.*(1 + n)))/4) .* dphi ...
             - (n .* (24 + n.*(24 + n*21))/8) .* sin(dphi)   .* cos(pphi) ...
                   + (15 * n.^2 .* (1 + n)/8) .* sin(2*dphi) .* cos(2*pphi) ...
                           - ((35 * n.^3)/24) .* sin(3*dphi) .* cos(3*pphi));

%-----------------------------------------------------------------------

function phiP = phicentral(y, ellipsoid, scalefactor, phi0, y0);
% Fixed point iteration for the latitude of the point (y,x0), which
% lies on the central meridian.

F0 = scalefactor;
a = ellipsoid(1);

tol = 1e-12 * a;

phiP = (y - y0)/(a * F0) + phi0;
yc = ycentral(phiP, ellipsoid, F0, phi0, y0);
nIterations = 0;
while any(abs(y - yc) > tol) && nIterations < 10
    phiP = (y - yc)/(a * F0) + phiP;
    yc = ycentral(phiP, ellipsoid, scalefactor, phi0, y0);
    nIterations = nIterations + 1;
end
