Index: /issm/trunk/src/m/utils/LatLong/ll2xy.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/ll2xy.m	(revision 2270)
+++ /issm/trunk/src/m/utils/LatLong/ll2xy.m	(revision 2271)
@@ -11,5 +11,5 @@
 re = 6378137.0; 	% WGS84
 e2 = 0.00669437999015;  % WGS84
-sn=-1.0; 	% because it's southern hemisphere
+sn=+1.0; 	% because it's southern hemisphere
 
 a=re;
Index: /issm/trunk/src/m/utils/LatLong/mapll.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/mapll.m	(revision 2270)
+++ /issm/trunk/src/m/utils/LatLong/mapll.m	(revision 2271)
@@ -1,3 +1,3 @@
-function  [x,y]=mapll(lat,lon,hem);
+function  [x,y]=mapll(lat,lon,hem,varargin);
 %MAPLL - convert latitude and longitude into x and y
 %
@@ -10,4 +10,37 @@
 %   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
@@ -17,17 +50,4 @@
 e2= 0.00669437999015;
 e=sqrt(e2);
-
-%Standard parallel - latitude with no distortion = -71.
-if hem==1,
-	%Northern hemisphere
-	sn=1;
-	slat=70;
-	xlam=45;
-else
-	%Southern Hemisphere
-	sn=-1.0;
-	slat=71;
-	xlam=0;
-end
 
 lat=sn*lat*pi/180;
@@ -48,2 +68,8 @@
 x= rho*sn.*sin((lon+xlam));
 y=-rho*sn.*cos((lon+xlam));
+
+end
+
+function mapllerrorusage(),
+	help mapll
+end
Index: /issm/trunk/src/m/utils/LatLong/mapxy.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/mapxy.m	(revision 2270)
+++ /issm/trunk/src/m/utils/LatLong/mapxy.m	(revision 2271)
@@ -1,3 +1,3 @@
-function [alat,alon]=mapxy(x,y,hem);
+function [alat,alon]=mapxy(x,y,hem,varargin);
 %MAPXY - compute latitude and longitude from x and y
 %
@@ -9,19 +9,40 @@
 %   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 ~(hem=='s' | hem=='n'),
+	error('mapxy 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
+
 re = 6378137;
 e2 = 0.00669437999015;
 
 e=sqrt(e2);
-
-%Standard parallel = latitude with no distortion
-slat=71;
-sn=-1.0;
-xlam=0;
-
-if hem ==1,
-   slat=70;
-	sn=1;
-	xlam=45;
-end
 
 slat=slat/180*pi;
@@ -54,2 +75,10 @@
    alon(indice)=alon(indice)-360;
 end
+
+
+end
+	
+
+function mapxyerrorusage()
+	help mapxy
+end
Index: /issm/trunk/src/m/utils/LatLong/stereomap.m
===================================================================
--- /issm/trunk/src/m/utils/LatLong/stereomap.m	(revision 2271)
+++ /issm/trunk/src/m/utils/LatLong/stereomap.m	(revision 2271)
@@ -0,0 +1,30 @@
+function new_data=stereomap(x_m,y_m,data,slat,xhem,new_x_m,new_y_m,new_slat,new_xhem,hemisphere);
+%STEREOMAP take a 2d map, on a stereographic projection defined by (slat, xhem), 
+%          and remap it to a new projection (on 2d axis (new_x_m,new_y_m) and new_slat,new_xhem,hemisphere)
+%
+%   usage: new_data=stereomap(x_m,y_m,data,slat,xhem,new_x_m,new_y_m,new_slat,new_xhem,hemisphere);
+%
+%   see also: mapll and mapxy
+
+M=length(new_y_m)-1;
+N=length(new_x_m)-1;
+
+%First, transform (new_x_m,new_y_m) into current projection: 
+x=zeros(M*N,1);
+y=zeros(M*N,1);
+for i=1:M,
+	x( (i-1)*N+1:i*N)=(new_x_m(1:end-1)+new_x_m(2:end))/2;
+	y( (i-1)*N+1:i*N)=(new_y_m(i)+new_y_m(i+1))/2*ones(N,1);
+end
+
+[lat,long]=mapxy(x,y,hemisphere,new_slat,new_xhem);
+[new_x,new_y]=mapll(lat,long,hemisphere,slat,xhem);
+
+%ok, we have the transformed coordinates in the current projection, go pick up the corresponding 
+%values.
+new_data_line=InterpFromGrid(x_m,y_m,data,new_x,new_y,NaN);
+
+new_data=zeros(M,N);
+for i=1:M,
+	new_data(i,:)=new_data_line( (i-1)*N+1:i*N);
+end