Index: ../trunk-jpl/src/m/coordsystems/ll2utm.m
===================================================================
--- ../trunk-jpl/src/m/coordsystems/ll2utm.m	(revision 0)
+++ ../trunk-jpl/src/m/coordsystems/ll2utm.m	(revision 18460)
@@ -0,0 +1,200 @@
+function [x,y,f]=ll2utm(lat,lon,datum)
+%LL2UTM Lat/Lon to UTM coordinates precise conversion.
+%   [X,Y]=LL2UTM2(LAT,LON) or LL2UTM([LAT,LON]) converts coordinates 
+%	LAT,LON (in degrees) to UTM X and Y (in meters). Default datum is WGS84.
+%
+%	LAT and LON can be scalars, vectors or matrix. Outputs X and Y will
+%	have the same size as inputs.
+%
+%	LL2UTM(LAT,LON,DATUM) uses specific DATUM for conversion. DATUM can be
+%	a string in the following list:
+%		'wgs84': World Geodetic System 1984 (default)
+%		'nad27': North American Datum 1927
+%		'clk66': Clarke 1866
+%		'nad83': North American Datum 1983
+%		'grs80': Geodetic Reference System 1980
+%		'int24': International 1924 / Hayford 1909
+%	or DATUM can be a 2-element vector [A,F] where A is semimajor axis (in
+%	meters)	and F is flattening of the user-defined ellipsoid.
+%
+%	[X,Y,ZONE]=LL2UTM(...) returns also computed UTM ZONE (negative value 
+%	stands for southern hemisphere points).
+%
+%
+%	XY = LL2UTM(...) or without any output argument returns a 2-column 
+%	matrix [X,Y].
+%
+%	Notice:
+%		- LL2UTM does not perform cross-datum conversion.
+%		- precision is near a millimeter.
+%
+%
+%	Reference:
+%		I.G.N., Projection cartographique Mercator Transverse: Algorithmes,
+%		   Notes Techniques NT/G 76, janvier 1995.
+%
+%   Author: Francois Beauducel, <beauducel@ipgp.fr>
+%   Created: 2003-12-02
+%   Updated: 2014-04-20
+
+
+%	Copyright (c) 2001-2014, François Beauducel, covered by BSD License.
+%	All rights reserved.
+%
+%	Redistribution and use in source and binary forms, with or without 
+%	modification, are permitted provided that the following conditions are 
+%	met:
+%
+%	   * Redistributions of source code must retain the above copyright 
+%	     notice, this list of conditions and the following disclaimer.
+%	   * Redistributions in binary form must reproduce the above copyright 
+%	     notice, this list of conditions and the following disclaimer in 
+%	     the documentation and/or other materials provided with the distribution
+%	                           
+%	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
+%	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
+%	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
+%	ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
+%	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+%	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+%	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
+%	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
+%	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
+%	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
+%	POSSIBILITY OF SUCH DAMAGE.
+
+% Available datums
+datums = [ ...
+	{ 'wgs84', 6378137.0, 298.257223563 };
+	{ 'nad83', 6378137.0, 298.257222101 };
+	{ 'grs80', 6378137.0, 298.257222101 };
+	{ 'nad27', 6378206.4, 294.978698214 };
+	{ 'int24', 6378388.0, 297.000000000 };
+	{ 'clk66', 6378206.4, 294.978698214 };
+];
+
+if nargin == 1
+	if size(lat,2) ~= 2
+		error('Single input argument must be a 2-column matrix [LAT,LON].')
+	end
+	lon = lat(:,2);
+	lat = lat(:,1);
+end
+	
+if nargin < 1
+	error('Not enough input arguments.')
+end
+
+if all([numel(lat),numel(lon)] > 1) && any(size(lat) ~= size(lon))
+	error('LAT and LON must be the same size or scalars.')
+end
+
+if nargin < 3
+	datum = 'wgs84';
+end
+
+if ischar(datum)
+	if ~any(strcmpi(datum,datums(:,1)))
+		error('Unkown DATUM name "%s"',datum);
+	end
+	k = find(strcmpi(datum,datums(:,1)));
+	A1 = datums{k,2};
+	F1 = datums{k,3};	
+else
+	if numel(datum) ~= 2
+		error('User defined DATUM must be a vector [A,F].');
+	end
+	A1 = datum(1);
+	F1 = datum(2);
+end
+
+% constants
+D0 = 180/pi;	% conversion rad to deg
+
+K0 = 0.9996;	% UTM scale factor
+X0 = 500000;	% UTM false East (m)
+
+B1 = A1*(1 - 1/F1);
+E1 = sqrt((A1*A1 - B1*B1)/(A1*A1));
+
+p1 = lat/D0;					% Phi = Latitude (rad)
+l1 = lon/D0;					% Lambda = Longitude (rad)
+F0 = round((l1*D0 + 183)/6);	% UTM zone
+P0 = 0/D0;
+L0 = (6*F0 - 183)/D0;			% UTM origin longitude (rad)
+Y0 = 1e7*(p1 < 0);				% UTM false northern (m)
+N = K0*A1;
+
+C = coef(E1,0);
+B = C(1)*P0 + C(2)*sin(2*P0) + C(3)*sin(4*P0) + C(4)*sin(6*P0) + C(5)*sin(8*P0);
+YS = Y0 - N*B;
+
+C = coef(E1,2);
+L = log(tan(pi/4 + p1/2).*(((1 - E1*sin(p1))./(1 + E1*sin(p1))).^(E1/2)));
+z = complex(atan(sinh(L)./cos(l1 - L0)),log(tan(pi/4 + asin(sin(l1 - L0)./cosh(L))/2)));
+Z = N.*C(1).*z + N.*(C(2)*sin(2*z) + C(3)*sin(4*z) + C(4)*sin(6*z) + C(5)*sin(8*z));
+xs = imag(Z) + X0;
+ys = real(Z) + YS;
+
+% outputs zone if needed: scalar value if unique, or vector/matrix of the
+% same size as x/y in case of crossed zones
+if nargout > 2
+	fu = unique(F0.*sign(lat));
+	if isscalar(fu)
+		f = fu;
+	else
+		f = F0;
+	end
+end
+
+if nargout < 2
+	x = [xs(:),ys(:)];
+else
+	x = xs;
+	y = ys;
+end
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function c = coef(e,m)
+%COEF Projection coefficients
+%	COEF(E,M) returns a vector of 5 coefficients from:
+%		E = first ellipsoid excentricity
+%		M = 0 for transverse mercator
+%		M = 1 for transverse mercator reverse coefficients
+%		M = 2 for merdian arc
+
+
+if nargin < 2
+	m = 0;
+end
+
+switch m
+	case 0
+	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
+           -105/4096, 0, -45/1024, 0,  -3/32, 0, -3/8, 0, 0;
+           525/16384, 0,  45/1024, 0, 15/256, 0,    0, 0, 0;
+          -175/12288, 0, -35/3072, 0,      0, 0,    0, 0, 0;
+          315/131072, 0,        0, 0,      0, 0,    0, 0, 0];
+	  
+	case 1
+	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
+             1/61440, 0,   7/2048, 0,   1/48, 0,  1/8, 0, 0;
+          559/368640, 0,   3/1280, 0,  1/768, 0,    0, 0, 0;
+          283/430080, 0, 17/30720, 0,      0, 0,    0, 0, 0;
+       4397/41287680, 0,        0, 0,      0, 0,    0, 0, 0];
+
+	case 2
+	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
+         -901/184320, 0,  -9/1024, 0,  -1/96, 0,  1/8, 0, 0;
+         -311/737280, 0,  17/5120, 0, 13/768, 0,    0, 0, 0;
+          899/430080, 0, 61/15360, 0,      0, 0,    0, 0, 0;
+      49561/41287680, 0,        0, 0,      0, 0,    0, 0, 0];
+   
+end
+c = zeros(size(c0,1),1);
+
+for i = 1:size(c0,1)
+    c(i) = polyval(c0(i,:),e);
+end
\ No newline at end of file
Index: ../trunk-jpl/src/m/coordsystems/utm2ll.m
===================================================================
--- ../trunk-jpl/src/m/coordsystems/utm2ll.m	(revision 18459)
+++ ../trunk-jpl/src/m/coordsystems/utm2ll.m	(revision 18460)
@@ -1,117 +1,190 @@
-function  [Lat,Lon] = utm2ll(xx,yy,utmzone)
-% -------------------------------------------------------------------------
-% [Lat,Lon] = utm2ll(x,y,utmzone)
+function [lat,lon]=utm2ll(x,y,f,datum,varargin)
+%UTM2LL UTM to Lat/Lon coordinates precise conversion.
+%   [LAT,LON]=UTM2LL(X,Y,ZONE) converts UTM coordinates X,Y (in meters)
+%	defined in the UTM ZONE (integer) to latitude LAT and longitude LON 
+%	(in degrees). Default datum is WGS84.
 %
-% Description: Function to convert vectors of UTM coordinates into Lat/Lon vectors (WGS84).
-% Some code has been extracted from UTMIP.m function by Gabriel Ruiz Martinez.
+%	X and Y can be scalars, vectors or matrix. Outputs LAT and LON will
+%	have the same size as inputs.
 %
-% Inputs:on)
-%    -3.485713    7.801235 -119.955246  -17.759537  -94.799019  121.640266
+%	For southern hemisphere points, use negative zone -ZONE.
 %
-% Example 2: If you need Lat/Lon coordinates in Degrees, Minutes and Seconds
-% [Lat, Lon]=utm2ll(x,y,utmzone);
-% LatDMS=dms2mat(deg2dms(Lat))
-%LatDMS =
-%    40.00         18.00         55.55
-%    46.00         17.00          2.01
-%    37.00         34.00         40.17
-%    28.00         38.00         44.33
-%    38.00         51.00         19.96
-%    25.00          3.00         42.41
-% LonDMS=dms2mat(deg2dms(Lon))
-%LonDMS =
-%    -3.00         29.00          8.61
-%     7.00         48.00          4.40
-%  -119.00         57.00         18.93
-%   -17.00         45.00         34.33
-%   -94.00         47.00         56.47
-%   121.00         38.00         24.96
+%	UTM2LL(X,Y,ZONE,DATUM) uses specific DATUM for conversion. DATUM can be
+%	a string in the following list:
+%		'wgs84': World Geodetic System 1984 (default)
+%		'nad27': North American Datum 1927
+%		'clk66': Clarke 1866
+%		'nad83': North American Datum 1983
+%		'grs80': Geodetic Reference System 1980
+%		'int24': International 1924 / Hayford 1909
+%	or DATUM can be a 2-element vector [A,F] where A is semimajor axis (in
+%	meters)	and F is flattening of the user-defined ellipsoid.
 %
-% Author:
-%   Rafael Palacios
-%   Universidad Pontificia Comillas
-%   Madrid, Spain
-% Version: Apr/06, Jun/06, Aug/06
-% Aug/06: corrected m-Lint warnings
-%-------------------------------------------------------------------------
+%	Notice:
+%		- UTM2LL does not perform cross-datum conversion.
+%		- precision is near a millimeter.
+%
+%
+%	Reference:
+%		I.G.N., Projection cartographique Mercator Transverse: Algorithmes,
+%		   Notes Techniques NT/G 76, janvier 1995.
+%
+%   Author: Francois Beauducel, <beauducel@ipgp.fr>
+%   Created: 2001-08-23
+%   Updated: 2014-04-20
 
-% Argument checking
+%	Copyright (c) 2001-2014, François Beauducel, covered by BSD License.
+%	All rights reserved.
 %
-error(nargchk(3, 3, nargin)); %3 arguments required
-n1=length(xx);
-n2=length(yy);
-n3=size(utmzone,1);
-if (n1~=n2 || n1~=n3)
-	error('x,y and utmzone vectors should have the same number or rows');
+%	Redistribution and use in source and binary forms, with or without 
+%	modification, are permitted provided that the following conditions are 
+%	met:
+%
+%	   * Redistributions of source code must retain the above copyright 
+%	     notice, this list of conditions and the following disclaimer.
+%	   * Redistributions in binary form must reproduce the above copyright 
+%	     notice, this list of conditions and the following disclaimer in 
+%	     the documentation and/or other materials provided with the distribution
+%	                           
+%	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
+%	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
+%	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
+%	ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
+%	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+%	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+%	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
+%	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
+%	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
+%	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
+%	POSSIBILITY OF SUCH DAMAGE.
+
+% Available datums
+datums = [ ...
+	{ 'wgs84', 6378137.0, 298.257223563 };
+	{ 'nad83', 6378137.0, 298.257222101 };
+	{ 'grs80', 6378137.0, 298.257222101 };
+	{ 'nad27', 6378206.4, 294.978698214 };
+	{ 'int24', 6378388.0, 297.000000000 };
+	{ 'clk66', 6378206.4, 294.978698214 };
+];
+
+if nargin < 3
+	error('Not enough input arguments.')
 end
-c=size(utmzone,2);
-if (c~=4)
-	error('utmzone should be a vector of strings like "30 T"');
+
+if all([numel(x),numel(y)] > 1) && any(size(x) ~= size(y))
+	error('X and Y must be the same size or scalars.')
 end
 
-% Memory pre-allocation
-%
-Lat=zeros(n1,1);
-Lon=zeros(n1,1);
+if ~isnumeric(f) || ~isscalar(f) || f ~= round(f)
+	error('ZONE must be a scalar integer.')
+end
 
-% Main Loop
-%
-for i=1:n1
-	if (utmzone(i,4)>'X' || utmzone(i,4)<'C')
-		fprintf('utm2ll: Warning utmzone should be a vector of strings like "30 T", not "30 t"\n');
+if nargin < 4
+	datum = 'wgs84';
+end
+
+if ischar(datum)
+	if ~any(strcmpi(datum,datums(:,1)))
+		error('Unkown DATUM name "%s"',datum);
 	end
-	if (utmzone(i,4)>'M')
-		hemis='N';   % Northern hemisphere
-	else
-		hemis='S';
+	k = find(strcmpi(datum,datums(:,1)));
+	A1 = datums{k,2};
+	F1 = datums{k,3};	
+else
+	if numel(datum) ~= 2
+		error('User defined DATUM must be a vector [A,F].');
 	end
+	A1 = datum(1);
+	F1 = datum(2);
+end
 
-	x=xx(i);
-	y=yy(i);
-	zone=str2double(utmzone(i,1:2));
+% constants
+D0 = 180/pi;	% conversion rad to deg
+maxiter = 100;	% maximum iteration for latitude computation
+eps = 1e-11;	% minimum residue for latitude computation
 
-	sa = 6378137.000000 ; sb = 6356752.314245;
+K0 = 0.9996;								% UTM scale factor
+X0 = 500000;								% UTM false East (m)
+Y0 = 1e7*(f < 0);							% UTM false North (m)
+P0 = 0;										% UTM origin latitude (rad)
+L0 = (6*abs(f) - 183)/D0;					% UTM origin longitude (rad)
+E1 = sqrt((A1^2 - (A1*(1 - 1/F1))^2)/A1^2);	% ellpsoid excentricity
+N = K0*A1;
 
-	%   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
+% computing parameters for Mercator Transverse projection
+C = coef(E1,0);
+YS = Y0 - N*(C(1)*P0 + C(2)*sin(2*P0) + C(3)*sin(4*P0) + C(4)*sin(6*P0) + C(5)*sin(8*P0));
 
-	X = x - 500000;
+C = coef(E1,1);
+zt = complex((y - YS)/N/C(1),(x - X0)/N/C(1));
+z = zt - C(2)*sin(2*zt) - C(3)*sin(4*zt) - C(4)*sin(6*zt) - C(5)*sin(8*zt);
+L = real(z);
+LS = imag(z);
 
-	if hemis == 'S' || hemis == 's'
-		Y = y - 10000000;
-	else
-		Y = y;
-	end
+l = L0 + atan(sinh(LS)./cos(L));
+p = asin(sin(L)./cosh(LS));
 
-	S = ( ( zone * 6 ) - 183 );
-	lat =  Y / ( 6366197.724 * 0.9996 );
-	v = ( c / ( ( 1 + ( e2cuadrada * ( cos(lat) ) ^ 2 ) ) ) ^ 0.5 ) * 0.9996;
-	a = X / v;
-	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 );
-	b = ( Y - Bm ) / v;
-	Epsi = ( ( e2cuadrada * a^ 2 ) / 2 ) * ( cos(lat) )^ 2;
-	Eps = a * ( 1 - ( Epsi / 3 ) );
-	nab = ( b * ( 1 - Epsi ) ) + lat;
-	senoheps = ( exp(Eps) - exp(-Eps) ) / 2;
-	Delt = atan(senoheps / (cos(nab) ) );
-	TaO = atan(cos(Delt) * tan(nab));
-	longitude = (Delt *(180 / pi ) ) + S;
-	latitude = ( lat + ( 1 + e2cuadrada* (cos(lat)^ 2) - ( 3 / 2 ) * e2cuadrada * sin(lat) * cos(lat) * ( TaO - lat ) ) * ( TaO - lat ) ) * ...
-		(180 / pi);
+L = log(tan(pi/4 + p/2));
 
-	Lat(i)=latitude;
-	Lon(i)=longitude;
+% calculates latitude from the isometric latitude
+p = 2*atan(exp(L)) - pi/2;
+p0 = NaN;
+n = 0;
+while any(isnan(p0) | abs(p - p0) > eps) & n < maxiter
+	p0 = p;
+	es = E1*sin(p0);
+	p = 2*atan(((1 + es)./(1 - es)).^(E1/2).*exp(L)) - pi/2;
+	n = n + 1;
+end
 
+if nargout < 2
+	lat = D0*[p(:),l(:)];
+else
+	lat = p*D0;
+	lon = l*D0;
 end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function c = coef(e,m)
+%COEF Projection coefficients
+%	COEF(E,M) returns a vector of 5 coefficients from:
+%		E = first ellipsoid excentricity
+%		M = 0 for transverse mercator
+%		M = 1 for transverse mercator reverse coefficients
+%		M = 2 for merdian arc
+
+
+if nargin < 2
+	m = 0;
+end
+
+switch m
+	case 0
+	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
+           -105/4096, 0, -45/1024, 0,  -3/32, 0, -3/8, 0, 0;
+           525/16384, 0,  45/1024, 0, 15/256, 0,    0, 0, 0;
+          -175/12288, 0, -35/3072, 0,      0, 0,    0, 0, 0;
+          315/131072, 0,        0, 0,      0, 0,    0, 0, 0];
+	  
+	case 1
+	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
+             1/61440, 0,   7/2048, 0,   1/48, 0,  1/8, 0, 0;
+          559/368640, 0,   3/1280, 0,  1/768, 0,    0, 0, 0;
+          283/430080, 0, 17/30720, 0,      0, 0,    0, 0, 0;
+       4397/41287680, 0,        0, 0,      0, 0,    0, 0, 0];
+
+	case 2
+	c0 = [-175/16384, 0,   -5/256, 0,  -3/64, 0, -1/4, 0, 1;
+         -901/184320, 0,  -9/1024, 0,  -1/96, 0,  1/8, 0, 0;
+         -311/737280, 0,  17/5120, 0, 13/768, 0,    0, 0, 0;
+          899/430080, 0, 61/15360, 0,      0, 0,    0, 0, 0;
+      49561/41287680, 0,        0, 0,      0, 0,    0, 0, 0];
+   
+end
+c = zeros(size(c0,1),1);
+
+for i = 1:size(c0,1)
+    c(i) = polyval(c0(i,:),e);
+end
\ No newline at end of file
