function [out1,out2,out3] = myutmcalc(in1,in2,direction,geoid,utm_zone_struct)

% Usage (radian angular units):
%
%	[easting,northing,utm_zone_struct] = myutmcalc(lon,lat,'forward',geoid)
%	[lon,lat] = myutmcalc(easting,northing,'inverse',geoid,utm_zone_struct)

%% ----------------------------------------------------------------------
%% unknown form of Transverse Mercator - Gauss-Kruger? (from Snyder's book)
%% ----------------------------------------------------------------------

% constants and parameters

k0 = .9996;							% scale on central meridian for UTM
a = geoid(1);						% semimajor axis in meters
e = geoid(2);						% eccentricity
e2 = e^2; e4 = e^4; e6 = e^6;		% eccentricity powers
ep2 = e2/(1-e2);					% second eccentricity parameter

if (in2 < 0), 
   falsenorthing = 10000000;
else
   falsenorthing = 0;
end;
   falseeasting = 500000;
      
% forward or inverse utm projection transform

switch direction

	case 'forward',   % forward transform
	
	  lon = in1;
	  lat = in2;
	
      out3 = myutmzone(median(lat),median(lon));
      phi  = (lat)*pi/180;
      dlam = (lon-out3.lon0)*pi/180;
     
      N = a./sqrt(1 - e2*(sin(phi)).^2);
      T = (tan(phi)).^2;
      C = ep2*(cos(phi)).^2;
      A = dlam.*cos(phi);
      
      Z1 = 1 - e2/4 - 3*e4/64 - 5*e6/256;
      Z2 = 3*e2/8 + 3*e4/32 + 45*e6/1024;
      Z3 = 15*e4/256 + 45*e6/1024;
      Z4 = 35*e6/3072;
		
      M = a*( Z1*phi - Z2*sin(2*phi) + Z3*sin(4*phi) - Z4*sin(6*phi) );
      % M is the true distance along the central meridian from the equator to lat
      M0 = 0; % Not used in this form of the calculation.
      % M0 is the true distance along the central meridian from the equator to the
      % central latitude for the particular UTM zone (ei, the origin latitude)
      %% lat_0 = out3.lat0*pi/180;
      %% M0 = a*( Z1*lat_0 - Z2*sin(2*lat_0) + Z3*sin(4*lat_0) - Z4*sin(6*lat_0) );
      
      Za = 1 - T + C;
      Zb = 5 - 18*T + T.^2 + 72*C - 58*ep2;
      Zc = 5 - T + 9*C + 4*C.^2;
      Zd = 61 - 58*T + T.^2 +600*C - 330*ep2;
      
      x = k0*N .* ( A + Za.*A.^3/6 + Zb.*A.^5/120 );
      y = k0 * ( M - M0 + N.*tan(phi).*( A.^2/2 + Zc.*A.^4/24 + Zd.*A.^6/720 ) );
      
      x = x + falseeasting;
      y = y + falsenorthing;
      out1 = x;
      out2 = y;
		
	case 'inverse',   % inverse transform

      if nargin < 5,
          disp('The UTM zone number and letter are required to covert from UTM to Lat/Lon');
          break; 
      end;
	  x = in1 - falseeasting;
	  y = in2 - falsenorthing;

	  e1 = (1-sqrt(1-e2))/(1+sqrt(1-e2));
	  M = y/k0;
	  mu = M / ( a*(1 - e2/4 - 3*e4/64 - 5*e6/256) );
		
	  Z1 = 3*e1/2 - 27*e1^3/32;
	  Z2 = 21*e1^2/16 - 55*e1^4/32;
	  Z3 = 151*e1^3/96;
	  Z4 = 1097*e1^4/512;
		
	  phi1 = mu + Z1*sin(2*mu) + Z2*sin(4*mu) + Z3*sin(6*mu) + Z4*sin(8*mu);

	  C1 = ep2*(cos(phi1)).^2;
	  T1 = (tan(phi1)).^2;
	  N1 = a ./ sqrt( 1 - e2*(sin(phi1)).^2 );
	  %R1 = a*(1-e2) ./ ( 1 - e2*(sin(phi1)).^2 ).*sqrt(( 1 - e2*(sin(phi1)).^2 ));
      R1 = a*(1-e2) ./ ( 1 - e2*(sin(phi1)).^2 ).^(3/2);
	  D = x ./ (N1*k0);

	  Za = 5 + 3*T1 + 10*C1 - 4*C1.^2 - 9*ep2;
	  Zb = 61 + 90*T1 + 298*C1 + 45*T1.^2 - 252*ep2 - 3*C1.^2;
	  Zc = 1 + 2*T1 + C1;
	  Zd = 5 - 2*C1 + 28*T1 - 3*C1.^2 + 8*ep2 + 24*T1.^2;
	  phi = phi1 - (N1.*tan(phi1)./R1).*(D.^2/2 - Za.*D.^4/24 + Zb.*D.^6/720);
	  dlam = (D - Zc.*D.^3./6 + Zd.*D.^5/120) ./ cos(phi1);

	  out1 = roundn(angledim(dlam,'radians','degrees')+utm_zone_struct.lon0,-7);
	  out2 = roundn(angledim(phi,'radians','degrees'),-7);
      
	otherwise,

		error('Unrecognized direction string')
		
end
