Index: sm/trunk/src/m/utils/LatLong/latlongfromxy.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/latlongfromxy.m	(revision 6340)
+++ 	(revision )
@@ -1,25 +1,0 @@
-function md=latlongfromxy(md,projection,hemisphere,merref,latref)
-% LATLONGFROMXY - create latitude and longitude of the mesh
-%
-%   projection can be 'stereographic' only
-%   hemisphere must be 'n' or 's'
-%   the last two arguments are the central meridien and the zero distorsion 
-%   latitude of the projection
-%   
-%   Usage:
-%      md=latlongfromxy(md,projection,hemisphere,merref,latref)
-%
-%   Example
-%      md=latlongfromxy(md,'stereographic','n',70,45)
-
-
-%check that the hemisphere is either 'n' or 's': 
-if ~(strcmpi(hemisphere,'n') | strcmpi(hemisphere,'s')),
-	error('latlongfromxy error message: hemisphere should be ''s'' or ''n''');
-end
-
-if ~strcmpi(projection,'stereographic'),
-	error('latlongfromxy error message: only stereographic projection implemented yet');
-end
-
-[md.lat,md.long]=mapxy(md.x,md.y,hemisphere,merref,latref);
Index: /issm/trunk/src/m/utils/LatLong/ll2xy.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/ll2xy.m	(revision 6340)
+++ /issm/trunk/src/m/utils/LatLong/ll2xy.m	(revision 6341)
@@ -1,30 +1,58 @@
-function [x,y]=ll2xy(lat,lon);
-%LL2XY - convert latitude and longitude coordinates to x and y
+function [x,y] = ll2xy(lat,lon,sgn)  
+%LL2XY - converts lat long to polar stereographic
+%
+%   Converts from geodetic latitude and longitude to Polar 
+%   Stereographic (X,Y) coordinates for the polar regions.
+%   Author: Michael P. Schodlok, December 2003 (map2ll)
 %
 %   Usage:
-%      [x,y]=ll2xy(lat,lon)
-%
-%   See also MAPLL, MAPXY
+%      [x,y] = map2ll(lat,lon,sgn)
+%      - sgn = Sign of latitude +1 : north latitude
+%                               -1 : south latitude
 
-lon=lon+360; % to have 0<lon<360
+if nargin ~= 3
+	help map2ll
+	error('bad usage');
+end
 
-re = 6378137.0; 	% WGS84
-e2 = 0.00669437999015;  % WGS84
-sn=-1.0; 	% because it's southern hemisphere
+% Conversion constant from degrees to radians
+cde  = 57.29577951;
+% Standard latitude for the SSM/I grid with no distorsion
+slat = 70.;
+% Radius of the earth in meters
+re   = 6378.273*10^3;
+% Eccentricity of the Hughes ellipsoid squared
+ex2   = .006693883;
+% Eccentricity of the Hughes ellipsoid
+ex    =  sqrt(ex2);
 
-a=re;
-e=sqrt(e2);
+if sgn == 1 
+	delta = 45.;
+else
+	delta = 0.0;
+end
 
-phi  = sn*lat*pi/180;
-lambda  = lon*pi/180; 
+latitude  = abs(lat) * pi/180.;
+longitude = (lon + delta) * pi/180.;
 
-qp = 1 - (1-e2)/2/e*log((1-e)/(1+e));
-%m=cos(phi)/sqrt(1-e2.*sin(phi)*sin(phi));
-q=(1-e2)*(sin(phi)./(1-e2*sin(phi).^2)-0.5/e*log((1-e*sin(phi))./(1+e*sin(phi))));
-rho = a*sqrt(qp-q);
-x = rho.*sin(lambda);
-y = - sn*rho.*cos(lambda);
+% compute X and Y in grid coordinates.
+T = tan(pi/4-latitude/2) ./ ((1-ex*sin(latitude))./(1+ex*sin(latitude))).^(ex/2);
 
+if (90 - slat) <  1.e-5 
+	rho = 2.*re*T/sqrt((1.+ex)^(1.+ex)*(1.-ex)^(1.-ex));
+else
+	sl  = slat*pi/180.;
+	tc  = tan(pi/4.-sl/2.)/((1.-ex*sin(sl))/(1.+ex*sin(sl)))^(ex/2.);
+	mc  = cos(sl)/sqrt(1.0-ex2*(sin(sl)^2));
+	rho = re*mc*T/tc;
+end
 
-%lon=lon-360; % to recover the initial longitude
+y = -rho .* sgn .* cos(sgn.*longitude);
+x =  rho .* sgn .* sin(sgn.*longitude);
 
+[cnt1,cnt2] = find(latitude >= pi / 2.);
+
+if cnt1
+	x(cnt1,1) = 0.0
+	y(cnt1,1) = 0.0
+end
Index: sm/trunk/src/m/utils/LatLong/mapll.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/mapll.m	(revision 6340)
+++ 	(revision )
@@ -1,82 +1,0 @@
-function  [x,y]=mapll(lat,lon,hem,varargin);
-%MAPLL - convert latitude and longitude into x and y
-%
-%   latitude and longitude are in degrees.
-%   Hemisphere must be 'n' for north and 's' for south
-%   One can specify the standard parallel and centre meridians.
-%   Default values are (70,45) for northern hemisphere and (71,0) for southerm hemisphere.
-%
-%   Usage:
-%      [x,y]=mapll(latitude,longitude,hemisphere)
-%      [x,y]=mapll(latitude,longitude,hemisphere,parallel,meridien)
-%
-%   Example:
-%      [x,y]=mapll(lat,lon,'s')
-%      [x,y]=mapll(lat,lon,'n',71,39)
-%
-%   See also MAPXY, LL2XY
-
-%check number of arguments.
-if ~((nargin==3) | (nargin==5)),
-	mapllerrorusage();
-end
-
-%hem: must be either 's' or 'n'
-if ~ischar(hem),
-	error('mapll error message: hemisphere argument should be ''n'' or ''s''');
-end
-
-%check hem is either 'n' or 's': 
-if ~(hem=='s' | hem=='n'),
-	error('mapll error message: hem should be ''s'' or ''n''');
-end
-
-%set sn: 
-if hem=='s',
-	slat=71;
-	sn=-1.0;
-	xlam=0;
-else
-	slat=70;
-	sn=1.0;
-	xlam=45;
-end
-
-%set defaults for standard parallels and centre meridians.
-if nargin==5,
-	slat=varargin{1};
-	xlam=varargin{2};
-end
-
-
-%some corrections
-lon=lon+360;       % to have 0<lon<360
-
-%WGS84
-re=6378137.0;
-e2= 0.00669437999015;
-e=sqrt(e2);
-
-lat=sn*lat*pi/180;
-lon=lon*pi/180;
-
-rlat=lat;
-slat=slat/180*pi;
-xlam=xlam/180*pi;
-
-t1=tan(pi/4.-rlat/2.)./((1.0-e*sin(rlat))./...
-(1.0+e*sin(rlat))).^(e/2.);
-
-t2=tan(pi/4.-slat/2.)./((1.0-e*sin(slat))./...
-(1.0+e*sin(slat))).^(e/2.);
-
-cm=cos(slat)./sqrt(1.0-e2*(sin(slat).^2));
-rho=re*cm.*t1./t2;
-x= rho*sn.*sin((lon+xlam));
-y=-rho*sn.*cos((lon+xlam));
-
-end
-
-function mapllerrorusage(),
-	help mapll
-end
Index: sm/trunk/src/m/utils/LatLong/mapxy.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/mapxy.m	(revision 6340)
+++ 	(revision )
@@ -1,91 +1,0 @@
-function [alat,alon]=mapxy(x,y,hem,varargin);
-%MAPXY - compute latitude and longitude from x and y
-%
-%   hemisphere must be 'n' for north and 's' for south
-%   One can specify the standard parallel and center meridien.
-%   Default values are (70,45) for northern hemisphere and (71,0) for southerm hemisphere.
-%
-%   Usage:
-%      [latitude,longitude]=mapxy(x,y,hemisphere)
-%      [latitude,longitude]=mapxy(x,y,hemisphere,parallel,meridien)
-%
-%   Example:
-%      [lat,long]=mapxy(x,y,'s')
-%      [lat,long]=mapxy(x,y,'n',71,39)
-%
-%   See also MAPLL, LL2XY
-
-%check number of arguments.
-if ~((nargin==3) | (nargin==5)),
-	mapxyerrorusage();
-end
-
-%hem: must be either 's' or 'n'
-if ~ischar(hem),
-	error('mapxy error message: hemisphere argument should be ''n'' or ''s''');
-end
-
-%check hem is either 'n' or 's': 
-if ~(strcmpi(hem,'n') | strcmpi(hem,'s')),
-	error('mapxy error message: hem should be ''s'' or ''n''');
-end
-
-%set sn: 
-if strcmpi(hem,'s'),
-	slat=71;
-	sn=-1.0;
-	xlam=0;
-else
-	slat=70;
-	sn=1.0;
-	xlam=45;
-end
-
-%set defaults for standard parallels and centre meridians.
-if nargin==5,
-	slat=varargin{1};
-	xlam=varargin{2};
-end
-
-re = 6378137;
-e2 = 0.00669437999015;
-
-e=sqrt(e2);
-
-slat=slat/180*pi;
-
-rho=sqrt(x.^2+y.^2);
-cm=cos(slat)./sqrt(1.0-e2*(sin(slat).^2));
-t=tan((pi/4)-(slat/2.))./((1.0-e*sin(slat))./(1.0+e*sin(slat))).^(e/2.);
-
-t=rho.*t./(re*cm);
-
-chi=(pi/2.)-2.*atan(t);
-
-alat=chi+((e2/2.)+(5.0*e2^2/24.)+(e2^3/12.))*sin(2*chi)+...
-((7.0*e2^2/48.)+(29.*e2^3/240.))*sin(4.0*chi)+...
-(7.0*e2^3/120.)*sin(6.0*chi);
-
-alat=(sn*alat/pi*180);
-xpr=sn*x;
-ypr=sn*y;
-alon=atan2(xpr,-ypr)/pi*180-sn*xlam;
-alon=sn*alon;
-
-indice=find(alon<0);
-if length(indice)>0,
-   alon(indice)=alon(indice)+360;
-end
-
-indice=find(alon>360);
-if length(indice)>0,
-   alon(indice)=alon(indice)-360;
-end
-
-
-end
-	
-
-function mapxyerrorusage()
-	help mapxy
-end
Index: /issm/trunk/src/m/utils/LatLong/utm2ll.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/utm2ll.m	(revision 6340)
+++ /issm/trunk/src/m/utils/LatLong/utm2ll.m	(revision 6341)
@@ -50,11 +50,8 @@
 end
 
-
-
 % Memory pre-allocation
 %
 Lat=zeros(n1,1);
 Lon=zeros(n1,1);
-
 
 % Main Loop
@@ -119,3 +116,2 @@
 
 end
-
Index: /issm/trunk/src/m/utils/LatLong/xy2ll.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/xy2ll.m	(revision 6341)
+++ /issm/trunk/src/m/utils/LatLong/xy2ll.m	(revision 6341)
@@ -0,0 +1,64 @@
+function [lat,lon] = xy2ll(x,y,sgn)
+%XY2LL - converts xy to lat long
+%
+%   Converts Polar  Stereographic (X,Y) coordinates for the polar regions to
+%   latitude and longitude Stereographic (X,Y) coordinates for the polar
+%   regions.
+%   Author: Michael P. Schodlok, December 2003 (map2xy.m)
+%
+%   Usage:
+%      [lat,lon] = map2xy(x,y,sgn);
+%      - sgn = Sign of latitude +1 : north latitude
+%                               -1 : south latitude
+
+if nargin ~= 3
+	help map2xy
+	error('bad usage');
+end
+
+% Conversion constant from degrees to radians
+cde  = 57.29577951;
+% Standard latitude for the SSM/I grid with no distorsion
+slat = 70.;
+% Radius of the earth in meters
+re   = 6378.273*10^3;
+% Eccentricity of the Hughes ellipsoid squared
+ex2   = .006693883;
+% Eccentricity of the Hughes ellipsoid
+ex    =  sqrt(ex2);
+
+if sgn == 1 
+	delta = 45.;
+else
+	delta = 0.0;
+end
+
+sl  = slat*pi/180.;
+rho = sqrt(x.^2 + y.^2);
+cm = cos(sl) / sqrt(1.0 - ex2 * (sin(sl)^2));
+T = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ex * sin(sl)) / (1.0 + ex * sin(sl)))^(ex / 2.0);
+
+if  abs(slat-90.) < 1.e-5
+	T = rho * sqrt((1. + ex)^(1. + ex) * (1. - ex)^(1. - ex)) / 2. / re;
+else
+	T = rho * T / (re * cm);
+end
+
+chi = (pi / 2.0) - 2.0 * atan(T);
+lat = chi + ((ex2 / 2.0) + (5.0 * ex2^2.0 / 24.0) + (ex2^3.0 / 12.0)) * ...
+	sin(2 * chi) + ((7.0 * ex2^2.0 / 48.0) + (29.0 * ex2^3 / 240.0)) * ...
+	sin(4.0 * chi) + (7.0 * ex2^3.0 / 120.0) * sin(6.0 * chi) ;
+
+lat = sgn * lat;
+lon = atan2(sgn * x,-sgn * y);
+lon = sgn * lon;
+
+[res1,res2] = find(rho <= 0.1);
+if res1
+	lat(res1,1) = 90. * sgn;
+	lon(res1,1) = 0.0;
+end
+
+lon = lon * 180. / pi;
+lat = lat * 180. / pi;
+lon = lon - delta; 
Index: sm/trunk/src/m/utils/LatLong/xyscale_south.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/xyscale_south.m	(revision 6340)
+++ 	(revision )
@@ -1,17 +1,0 @@
-function res=xyscale_south(lat);
-%XYSCALE_SOUTH - calculate scaling factor
-%
-%   latitude is <0 in degrees.
-%
-%   Usage:
-%      res=xyscale_south(latatitude)
-lat=-lat/180*pi;
-
-m71_t71=1.9390295644;
-e2=.00669437999015;
-e=sqrt(e2);
-
-m=cos(lat)./sqrt(1-e2*sin(lat).^2);
-t=tan(pi/4-lat/2)./((1-e*sin(lat))./(1+e*sin(lat))).^(e/2);
-k=m71_t71*t./m;
-res=(1./k./k);