Index: /issm/trunk-jpl/externalpackages/dem/dem.m
===================================================================
--- /issm/trunk-jpl/externalpackages/dem/dem.m	(revision 23662)
+++ /issm/trunk-jpl/externalpackages/dem/dem.m	(revision 23663)
@@ -1,75 +1,225 @@
-function [h,I]=dem(x,y,z,varargin)
+function varargout=dem(x,y,z,varargin)
 %DEM Shaded relief image plot
+%
 %	DEM(X,Y,Z) plots the Digital Elevation Model defined by X and Y 
 %	coordinate vectors and elevation matrix Z, as a lighted image using
 %	specific "landcolor" and "seacolor" colormaps. DEM uses IMAGESC 
 %	function which is much faster than SURFL when dealing with large 
-%	high-resolution DEM.
-%
-%	DEM(X,Y,Z,OPT) specifies options with OPT = [A,C,SCUT,ZMIN,ZMAX,ZCUT],
-%	where sorted optional scalars are:
-%	  A = azimuth light (in degrees relative to North). Default is A = -45
-%	      for a natural northwestern illumination.
-%	  C = controls contrast, as the exponent of the gradient value. Default
-%	      is C = 1 for linear contrast; use C = 0 to remove lighting, 
-%	      C = 0.5 for moderate lighting, C = 2 or more for strong contrast.
-%	  SCUT = controls lighting scale saturation with a median-style filter
-%	      in % of elements. Default is SCUT = 0.2 (0.2% maximum gradient 
-%	      values is ignored). Use SCUT = 0 for full scale gradient.
-%	  ZMIN,ZMAX = fixes min and max elevation values for colormap. Use NaN
-%	      to keep real min and/or max data values.
-%	  ZCUT = median-style filter to cut extremes values (in % of elements).
-%	      Default is ZCUT = 0.5 to ignore the 0.5% of most min/max 
-%	      elevation values. Use ZCUT = 0 for full scale.
-%	Use OPT=[] to keep default values and define following arguments.
-%
-%	DEM(...,CMAP) uses CMAP colormap instead of default (landcolor, if 
-%	exists or jet). Note that defining CMAP disables the default SEACOLOR
-%	colormap attribution for Z<=0 elevations.
-%
-%	DEM(...,NOVALUE) defines the values that will be replaced by NaN. This 
-%	might be mandatory for DEM that use a value like -99999 or -32768 if 
-%	you don't want a flat image...
-%
-%	DEM(...,SEACOLOR) sets the colormap used for zero and negative values.
-%	Default is seacolor (if exists) or single color [0.7,0.9,1] (a light 
-%	cyan) to simulate sea color. Use [] to apply colormap CMAP on the full
-%	elevation scale.
-%
-%	DEM(...,'interp') interpolates linearly NaN values (fills the gaps).
-%
-%	DEM(...,'lake') detects automaticaly flat areas different from sea 
-%	level (non-zero elevations) and draws them as lake surfaces.
-%
-%	DEM(...,'dec') plots classic basemap-style axis, considering 
-%	coordinates as cartesian (decimal).
-%
-%	DEM(...,'dms') plots geographic basemap-style axis in deg/min/sec,
-%	considering coordinates X as longitude and Y as latitude. Axis aspect 
-%	ratio will be adjusted to approximatively preserve distances (this is  
-%	not a real projection!).
-%
-%	DEM(...,'scale') adds a legend to the right of graph, with elevation 
-%	scale (colormap) and a distance scale if 'dms' option is used.
-%
-%	[H,I]=DEM(...) returns graphic handle H and illuminated image as I, an
-%	MxNx3 matrix (if Z is MxN).
-%
-%	Informations:
-%	 - For optimization purpose, DEM will automatically decimate data to
-%	   limit to a total of 1500x1500 pixels images. To avoid it, use option
-%	   DEM(...,'nodecim') or DEM(...,'decim',N) where N is an integer, but 
-%	   be aware that large grids may require computer ressources and induce  
-%	   disk swap or memory errors.
-%	 - Colormaps are Mx3 RGB matrix so it is easy to modify saturation 
-%	   (CMAP.^N), set darker (CMAP/N), lighter ((N - 1 + CMAP)/N), inverse
-%	    it (flipud(CMAP)), etc...
-%	 - To get free worldwide topographic data (SRTM), see READHGT function.
+%	high-resolution DEM. It produces also high-quality and moderate-size 
+%	Postscript image adapted for publication.
+%
+%	DEM(X,Y,Z,'Param1',Value1,'Param2',Value2,...) specifies options or
+%	parameter/value couple (case insensitive):
+%
+%	[H,I] = DEM(...); returns graphic handle H and optional illuminated 
+%	image as I, a MxNx3 matrix (if Z is MxN and DECIM is 1).
+%
+%	I = DEM(...,'noplot') returns a structure I containing fields x, y, z,
+%	and illuminated image .rgb without producing a graph on current figure.
+%
+%
+%	--- Lighting options ---
+%
+%	'Azimuth',A
+%		Light azimuth in degrees clockwise relative to North. Default is
+%		A = -45 for	a natural northwestern illumination.
+%
+%	'Contrast',C
+%		Light contrast, as the exponent of the gradient value:
+%			C = 1 for linear contrast (default),
+%			C = 0 to remove lighting,
+%			C = 0.5 for moderate lighting,
+%			C = 2 or more for strong contrast.
+%
+%	'LCut',LC
+%		Lighting scale saturation cut with a median-style filter in % of 
+%	    elements, such as LC% of maximum gradient values are ignored:
+%			LC = 0.2 is default, 
+%			LC = 0 for full scale gradient.
+%
+%	'km'
+%		Stands that X and Y coordinates are in km instead of m (default).
+%		This allows correct lighting. Ignored if LATLON option is used.
+%
+%
+%	--- Elevation colorscale options ---
+%
+%	'ZLim',[ZMIN,ZMAX]
+%		Fixes min and max elevation values for colormap. Use NaN to keep 
+%		real min and/or max data values.
+%
+%	'ZCut',ZC
+%		Median-style filter to cut extremes values of Z (in % of elements),
+%		such that ZC% of most min/max elevation values are ignored in the
+%		colormap application:
+%			ZC = 0.5 is default, 
+%			ZC = 0 for full scale.
+%
+%
+%	--- "No Value" elevation options ---
+%
+%	'NoValue',NOVALUE
+%		Defines the values that will be replaced by NaN. Note that values 
+%		equal to minimum of Z class are automatically detected as NaN 
+%		(e.g., -32768 for int16 class).
+%
+%	'NaNColor',[R,G,B]
+%		Sets the RGB color for NaN/NoValue pixels (default is a dark gray).
+%		Note that your must specify a valid 3-scalar vector (between 0 and
+%		1);	color characters like 'w' or 'k' are not allowed, use [1,1,1]
+%		or [0,0,0] instead.
+%
+%	'Interp'
+%		Interpolates linearly all NaN values (fills the gaps using linear 
+%		triangulation), using an optimized algorithm.
+%
+%
+%	--- Colormap options ---
+%
+%	'LandColor',LMAP
+%		Uses LMAP colormap instead of default (landcolor, if exists or 
+%		jet) for Z > 0 elevations.
+%
+%	'SeaColor',SMAP
+%		Sets the colormap used for Z <= 0 elevations. Default is seacolor 
+%		(if exists) or single color [0.7,0.9,1] (a light cyan) to simulate
+%		sea color.
+%
+%	'ColorMap',CMAP
+%		Uses CMAP colormap for full range of elevations, instead of default 
+%		land/sea. This option overwrites LANDCOLOR/SEACOLOR options.
+%
+%	'Lake'
+%		Detects automaticaly flat areas different from sea level (non-zero 
+%		elevations) and colors them as lake surfaces.
+%
+%	'Watermark',N
+%		Makes the whole image lighter by a factor of N.
+%
+%
+%	--- Basemap and scale options ---
+%
+%	'Legend'
+%		Adds legends to the right of graph: elevation scale (colorbar)
+%		and a distance scale (in km).
+%
+%	'Cartesian'
+%		Plots classic basemap-style axis, considering coordinates X and Y 
+%		as cartesian in meters. Use parameter "km' for X/Y in km.
+%
+%	'LatLon'
+%		Plots geographic basemap-style axis in deg/min/sec, considering 
+%		coordinates X as longitude and Y as latitude. Axis aspect ratio 
+%		will be adjusted to approximatively preserve distances (this is  
+%		not a real projection!). This overwrites ZRatio option.
+%
+%	'AxisEqual', 'auto' (default) | 'manual' | 'off'
+%		When 'Cartesian' or 'LatLon' option is used, automatic axes scaling
+%		is applied to respect data aspect ratio. Default mode is 'auto' and
+%		uses AXIS EQUAL and DASPECT functions. The 'manual' mode modifies
+%		axes width or height with respect to the paper size in order to
+%		produce correct data scaling at print (but not necessarily at 
+%		screen). The 'off' mode disables any scaling.
+%
+%	Additionnal options for basemap CARTESIAN or LATLON:
+%
+%	'BorderWidth',BW
+%		Border width of the basemap axis, in % of axis height. Default is
+%		BW = 1%.
+%
+%	'XTick',DX
+%	'YTick',DY
+%		X and Y Tick length (same unit as X and Y). Default is automatic.
+%		Tick labels are every 2 ticks.
+%
+%	'FontSize',FS
+%		Font size for X and Y tick labels. Default is FS = 10.
+%
+%	'FontBold'
+%		Font weight bold for tick labels.
+%
+%	'Position',P
+%		Position of the tick labels: 'southwest' (default), 'southeast',
+%		'northwest','northeast'
+%
+%
+%	--- Decimation options ---
+%
+%	For optimization purpose, DEM will automatically decimate data to limit
+%	to a total of 1500x1500 pixels images. To avoid this, use following
+%	options, but be aware that large grids may require huge computer 
+%	ressources or induce disk swap or memory errors.
+%
+%	'Decim',N
+%		Decimates matrix Z at 1/N times of the original sampling.
+%
+%	'NoDecim'
+%		Forces full resolution of Z, no decimation.
+%
+%
+%
+%	--- Informations ---
+%
+%	Colormaps are Mx3 RGB matrix so it is easy to modify saturation 
+%	(CMAP.^N), set darker (CMAP/N), lighter (1 - 1/N + CMAP/N), inverse
+%	it (flipud(CMAP)), etc...
+%
+%	To get free worldwide topographic data (SRTM), see READHGT function.
+%
+%	For backward compatibility, the former syntax is still accepted:
+%	DEM(X,Y,Z,OPT,CMAP,NOVALUE,SEACOLOR) where OPT = [A,C,LC,ZMIN,ZMAX,ZC],
+%	also option aliases DEC, DMS and SCALE, but there is no argument 
+%	checking. Please prefer the param/value syntax.
 %
 %	Author: François Beauducel <beauducel@ipgp.fr>
-%	Created: 2007-05-17
-%	Updated: 2013-01-05
-
-%	Copyright (c) 2013, François Beauducel, covered by BSD License.
+%	Created: 2007-05-17 in Guadeloupe, French West Indies
+%	Updated: 2016-01-31
+
+%	History:
+%	[2016-04-19] v2.3
+%		- major update (thanks to mas Wiwit)
+%	[2016-01-31] v2.2
+%		- adds option 'Position' for tick labels
+%	[2015-08-22] v2.1
+%		- minor fix (former versions of Matlab compatibility)
+%	[2015-08-19] v2.0
+%		- image is now 100% true color (including the legend colorbar), 
+%	      thus completely independent from the figure colormap
+%	[2014-10-14]
+%		- 'decim' option allows oversampling (negative value)
+%	[2014-06-06]
+%		- improves backward compatibility (adds strjoin subfunction)
+%	[2014-03-18]
+%		- adds new axisequal option
+%	[2013-03-11]
+%		- new options: 'km', 'watermark', 'fontsize', 'bordersize'
+%		- improve legend colorbar
+%		- all options now passed as param/value
+%	[2013-01-14]
+%		- improved light rendering (using surface normals instead of gradient)
+%		- improved 'lake' detection algorithm
+%		- new 'nancolor' option to set NaN color
+%		- adds a length scale with 'dec' option
+%		- minor code improvements
+%	[2013-01-07]
+%		- adds 'interp' option (fill the gaps)
+%		- adds 'seacolor' colormap for negative elevations (bathymetry)
+%	[2013-01-02]
+%		- adds a 'lake' option
+%		- minor bug correction
+%	[2012-09-26]
+%		- now accepts row/column vectors for X and/or Y.
+%	[2012-05-29]
+%		- adds basemap-style axis in decimal or lat/lon modes
+%		- adds elevation and distance scales
+%	[2012-05-18]
+%		- new landcolor.m colormap function
+%		- new arguments to control colormap scaling
+%		- median-style filters for light and colormap
+%	[2012-04-26]
+%		- Optimizations: adds a decimation for large DEM grids.
+
+%
+%	Copyright (c) 2016, François Beauducel, covered by BSD License.
 %	All rights reserved.
 %
@@ -100,80 +250,233 @@
 end
 
-% default OPT arguments
-a = -45;
-c = 1;
-scut = 0.2;
-zmin = NaN;
-zmax = NaN;
-zcut = 0.5;
-grey = 0.2*[1,1,1];
-csea = [];
-fs = 10;	% tick label fontsize
-
-dec = 0;
-dms = 0;
-scale = 0;
-inter = 0;
-lake = 0;
-decimflag = 0;
-decim = 0;
-
-novalue_color = [0,0,0];
-
-if ~isnumeric(x) | ~isnumeric(y) | ~isnumeric(z)
+holdon = ishold;
+
+degkm = 6378*pi/180; % one latitude degree in km
+sea_color = [.7,.9,1]; % default sea color (light cyan)
+grey = 0.2*[1,1,1]; % a dark gray
+
+
+% -------------------------------------------------------------------------
+% --- Manage input arguments
+
+% number of arguments param/value
+nargs = 0;
+
+if ~isnumeric(x) || ~isnumeric(y) || ~isnumeric(z)
 	error('X,Y and Z must be numeric.')
 end
 
-if all(size(x) ~= 1) | all(size(y) ~= 1)
+if all(size(x) ~= 1) || all(size(y) ~= 1)
 	error('X and Y must be vectors, not matrix.')
 end
 
-if length(x) ~= size(z,2) | length(y) ~= size(z,1)
+if length(x) ~= size(z,2) || length(y) ~= size(z,1)
 	error('If Z has a size of [M,N], X must have a length of N, and Y a length of M.')
 end
 
-if nargin > 3
-	dec = any(strcmp(varargin,'dec'));
-	dms = any(strcmp(varargin,'dms'));
-	if dms & any(abs(y) > 91)
-		error('With DMS option Y must be in valid latitudes interval (decimal degrees).')
-	end
-	scale = any(strcmp(varargin,'scale'));
-	inter = any(strcmp(varargin,'interp'));
-	lake = any(strcmp(varargin,'lake'));
-	if any(strcmp(varargin,'nodecim'))
-		decim = 1;
-		decimflag = 1;
-	end
-	kdecim = find(strcmp(varargin,'decim'));
-	if ~isempty(kdecim)
-		decimflag = 1;
-		if (kdecim + 1) <= nargin & isnumeric(varargin{kdecim+1})
-			decim = round(varargin{kdecim+1});
-			decimflag = 2;
-		end
-	end
-end
-nargs = decimflag + dec + dms + scale + lake + inter;
-
-if (nargin - nargs) > 3
+% OPTIONS and PARAM/VALUE arguments
+			
+% AZIMUTH param/value
+[s,az] = checkparam(varargin,'azimuth',@isscalar);
+nargs = nargs + 2;
+if s==0
+	az = -45; % default
+end
+
+% ELEVATION param/value
+[s,el] = checkparam(varargin,'elevation',@isscalar);
+nargs = nargs + 2;
+if s==0
+	el = 0; % default
+end
+
+% CONTRAST param/value
+[s,ct] = checkparam(varargin,'contrast',@isscalar);
+nargs = nargs + 2;
+if s
+	ct = abs(ct);
+else
+	ct = 1; % default
+end
+
+% LCUT param/value
+[s,lcut] = checkparam(varargin,'lcut',@isperc);
+nargs = nargs + 2;
+if s==0
+	lcut = .2; % default
+end
+
+% NOVALUE param/value
+[s,novalue] = checkparam(varargin,'novalue',@isscalar);
+nargs = nargs + 2;
+if s==0
+	% default: min value for integer class / NaN for float
+	S = whos('z');
+	if strfind(S.class,'int')
+		novalue = intmin(S.class);
+	else
+		novalue = NaN;
+	end
+end
+
+% NANCOLOR param/value
+[s,novalue_color] = checkparam(varargin,'nancolor',@isrgb);
+nargs = nargs + 2;
+if s==0
+	novalue_color = grey; % default
+end
+
+% LANDCOLOR param/value
+[s,cland] = checkparam(varargin,'landcolor',@isrgb);
+nargs = nargs + 2;
+if s==0
+	% default: landcolor or jet
+	if exist('landcolor','file')
+		cland = landcolor.^1.3;
+	else
+		cland = jet(256);
+	end
+end
+
+% SEACOLOR param/value
+[s,csea] = checkparam(varargin,'seacolor',@isrgb);
+nargs = nargs + 2;
+if s==0
+	% default: seacolor or single color
+	if exist('seacolor','file')
+		csea = seacolor;
+	else
+		csea = sea_color;
+	end
+end
+
+% COLORMAP param/value
+[s,cmap] = checkparam(varargin,'colormap',@isrgb);
+nargs = nargs + 2;
+if s
+	cland = [];
+	csea = [];
+else
+	% default
+	cmap = cland;
+end
+
+% ZLIM param/value
+[s,zmm] = checkparam(varargin,'zlim',@isvec);
+nargs = nargs + 2;
+if s
+	zmin = min(zmm);
+	zmax = max(zmm);
+else
+	zmin = NaN; % default
+	zmax = NaN; % default
+end
+
+% ZCUT param/value
+[s,zcut] = checkparam(varargin,'zcut',@isperc);
+nargs = nargs + 2;
+if s==0
+	zcut = .5; % default
+end
+
+% ZRATIO param/value
+[s,zratio] = checkparam(varargin,'zratio',@isscalar);
+nargs = nargs + 2;
+if s==0
+	zratio = 1; % default
+end
+
+% WATERMARK param/value
+[s,wmark] = checkparam(varargin,'watermark',@isscalar);
+nargs = nargs + 2;
+if s
+	wmark = abs(wmark);
+else
+	wmark = 0; % default
+end
+
+% DECIM param/value and NODECIM option
+[s,decim] = checkparam(varargin,'decim',@isscalar);
+if s
+	decim = round(decim);
+	nargs = nargs + 2;
+else
+	decim = any(strcmpi(varargin,'nodecim')); % default
+	nargs = nargs + 1;
+end
+
+% FONTSIZE param/value
+[s,fs] = checkparam(varargin,'fontsize',@isscalar);
+nargs = nargs + 2;
+if s==0
+	fs = 10; % default
+end
+
+% BORDERWIDTH param/value
+[s,bw] = checkparam(varargin,'borderwidth',@isperc);
+nargs = nargs + 2;
+if s==0
+	bw = 1; % default
+end
+
+% XTICK param/value
+[s,ddx] = checkparam(varargin,'xtick',@isscalar);
+nargs = nargs + 2;
+if s==0
+	ddx = 0; % default (automatic)
+end
+
+% YTICK param/value
+[s,ddy] = checkparam(varargin,'ytick',@isscalar);
+nargs = nargs + 2;
+if s==0
+	ddy = 0; % default (automatic)
+end
+
+% POSITION param/value
+[s,tpos] = checkparam(varargin,'position',@ischar,{'southwest','southeast','northwest','northeast'});
+nargs = nargs + 2;
+if s==0
+	tpos = 'southwest'; % default
+end
+
+% AXISEQUAL param/value
+[s,axeq] = checkparam(varargin,'axisequal',@ischar,{'auto','manual','off'});
+nargs = nargs + 2;
+if s==0 || ~any(strcmpi(axeq,{'manual','off'}))
+	axeq = 'auto'; % default (automatic)
+end
+
+% CROP param/value
+[s,crop] = checkparam(varargin,'crop',@isvec,4);
+nargs = nargs + 2;
+
+% options without argument value
+km = any(strcmpi(varargin,'km'));
+dec = any(strcmpi(varargin,'cartesian') | strcmpi(varargin,'dec'));
+dms = any(strcmpi(varargin,'latlon') | strcmpi(varargin,'dms'));
+scale = any(strcmpi(varargin,'legend') | strcmpi(varargin,'scale'));
+inter = any(strcmpi(varargin,'interp'));
+lake = any(strcmpi(varargin,'lake'));
+fbold = any(strcmpi(varargin,'fontbold'));
+noplot = any(strcmpi(varargin,'noplot'));
+
+
+% for backward compatibility (former syntax)...
+nargs = nargs + dec + dms + scale + lake + inter + km + fbold + noplot;
+
+if (nargin - nargs) > 3 && ~isempty(varargin{1})
 	opt = varargin{1};
 	if ~isnumeric(opt)
-		error('OPT = [A,C,S,ZMIN,ZMAX] argument must be numeric.');
+		error('OPT = [A,C,S,ZMIN,ZMAX,ZCUT] argument must be numeric.');
 	end
 	if ~isempty(opt)
-		a = opt(1);
+		az = opt(1);
 	end
 	if length(opt) > 1
-		c = opt(2);
-		if c < 0
-			error('C argument must be positive.');
-		end
+		ct = opt(2);
 	end
 	if length(opt) > 2
-		scut = opt(3);
-		if scut < 0 | scut >= 100
-			error('SCUT argument must be a positive percentage.');
-		end
+		lcut = opt(3);
 	end
 	if length(opt) > 4
@@ -183,54 +486,44 @@
 	if length(opt) > 5
 		zcut = opt(6);
-		if zcut < 0 | zcut >= 100
-			error('ZCUT argument must be a positive percentage.');
-		end
-	end
-end
-
-if (nargin - nargs) < 5
-	cmap = [];
-else
+	end
+end
+
+if (nargin - nargs) > 4 && ~isempty(varargin{2})
 	cmap = varargin{2};
-	if ~isnumeric(cmap) | (~isempty(cmap) & (size(cmap,2) ~= 3 | min(cmap(:)) < 0 | max(cmap(:)) > 1))
-		error('CMAP must be a valid colormap (3-column [R,G,B] matrix with 0.0 to 1.0 values).')
-	end
-end
-
-if (nargin - nargs) < 6
-	novalue = NaN;
-else
+	csea = [];
+end
+
+if (nargin - nargs) > 5 && ~isempty(varargin{3})
 	novalue = varargin{3};
-	if ~isnumeric(novalue) | numel(novalue) > 1
-		error('NOVALUE must be scalar.')
-	end
-end
-
-if (nargin - nargs) < 7
-	if isempty(cmap)
-		if exist('seacolor','file')
-			csea = seacolor(256);
-		else
-			csea = [.7,.9,1];
-		end
-	end
-else
+end
+
+if (nargin - nargs) > 6 && ~isempty(varargin{4})
 	csea = varargin{4};
-	if ~isnumeric(csea)
-		error('Unknown option')
-	elseif (~isempty(csea) & (size(csea,2) ~= 3 | min(csea) < 0 | max(csea) > 1))
-		error('SEACOLOR must be a valid [R,G,B] vector with 0.0 to 1.0 values).')
-	end
-end
-
-if isempty(cmap)
-	if exist('landcolor','file')
-		cmap = landcolor.^1.3;
-	else
-		cmap = jet(256);
-	end
-end
-
-
+end
+
+
+% further test of input arguments
+if dms && any(abs(y) > 91)
+	error('With LATLON option Y must be in valid latitudes interval (decimal degrees).')
+end
+
+if km
+	zratio = 1000;
+end
+
+
+% -------------------------------------------------------------------------
+% --- Pre-process DEM data
+
+% crops data if needed
+if numel(crop)==4
+	fprintf('DEM: crops original data from [%g,%g,%g,%g] to [%g,%g,%g,%g]...\n', ...
+		min(x(:)),max(x(:)),min(y(:)),max(y(:)),crop);
+	kx = find(x >= crop(1) & x <= crop(2));
+	ky = find(y >= crop(3) & y <= crop(4));
+	x = x(kx);
+	y = y(ky);
+	z = z(ky,kx);
+end
 
 % decimates data to avoid disk swap/out of memory...
@@ -245,5 +538,5 @@
 	y = y(1:n:end);
 	z = z(1:n:end,1:n:end);
-	fprintf('DEM: on the plot data has been decimated by a factor of %d...\n',n);
+	fprintf('DEM: data has been decimated by a factor of %d...\n',n);
 end
 
@@ -251,74 +544,98 @@
 z(z==novalue) = NaN;
 
+if isempty(csea)
+	k = (z~=0 & ~isnan(z));
+else
+	k = ~isnan(z);
+end
+
+if isnan(zmin)
+	zmin = nmedian(z(k),zcut/100);
+end
+if isnan(zmax)
+	zmax = nmedian(z(k),1 - zcut/100);
+end
+dz = zmax - zmin;
+
+if decim & n < 0
+	xi = linspace(x(1),x(end),-n*length(x));
+	yi = linspace(y(1),y(end),-n*length(y))';
+	[xx,yy] = meshgrid(xi,yi);
+	z = interp2(x,y,z,xx,yy,'*cubic');
+	x = xi;
+	y = yi;
+	fprintf('DEM: data has been oversampled by a factor of %d...\n',-n);
+end
+
 if inter
 	z = fillgap(x,y,z);
 end
 
-if isempty(csea)
-	k = (z~=0 & ~isnan(z));
-else
-	k = ~isnan(z);
-end
-
-if isnan(zmin)
-	zmin = nmedian(z(k),zcut/100);
-end
-if isnan(zmax)
-	zmax = nmedian(z(k),1 - zcut/100);
-end
-dz = zmax - zmin;
+% -------------------------------------------------------------------------
+% --- Process lighting
 
 if dz > 0
 	% builds the colormap: concatenates seacolor and landcolor around 0
-	if ~isempty(csea) & zmin < 0 & zmax > 0
+	if ~isempty(csea)
 		l = size(csea,1);
-		r = size(cmap,1)*abs(zmin)/zmax/l;
-		cmap = cat(1,interp1(1:l,csea,linspace(1,l,round(l*r)),'*linear'),cmap);
+		if zmin < 0 && zmax > 0
+			r = size(cland,1)*abs(zmin)/zmax/l;
+			cmap = cat(1,interp1(1:l,csea,linspace(1,l,ceil(l*r)),'*linear'),cland);
+		elseif zmax <=0
+			cmap = csea;
+		end
 	end
 	
 	% normalisation of Z using CMAP and convertion to RGB
-	I = ind2rgb(uint16((z - zmin)*(length(cmap)/dz)),cmap);
+	I = ind2rgb(uint16(round((z - zmin)*(size(cmap,1) - 1)/dz) + 1),cmap);
 	
-	if c > 0
+	if ct > 0
 		% computes lighting from elevation gradient
-		[fx,fy] = gradient(z,x,y);
-		%fx = filter([1,0,0,0,-1],1,rf(z'))';
-		%fy = filter([1,0,0,0,-1],1,rf(z));
-		fxy = -fx*sind(a) - fy*cosd(a);
-		clear fx fy	% free some memory...
+		%[fx,fy] = gradient(z,x,y);
+		if dms
+			ryz = degkm*1000;
+			rxz = degkm*1000*cosd(mean(y));
+		else
+			rxz = zratio;
+			ryz = zratio;
+		end
+		[xx,yy] = meshgrid(x*rxz,y*ryz);
+		[fx,fy,fz] = surfnorm(xx,yy,z);
+		[ux,uy,uz] = sph2cart((90-az)*pi/180,el*pi/180,1);
+		fxy = fx*ux + fy*uy + fz*uz;
+		clear xx yy fx fy fz	% free some memory...
 		
-		% lake option: zero gradient
-		if lake
-			dx = diff(z,1,2);
-			dy = diff(z,1,1);
-			u1 = ones(size(z,1),1);
-			u2 = ones(1,size(z,2));
-			z(cat(1,u2,dy)==0 & cat(1,dy,u2)==0 & cat(2,dx,u1)==0 & cat(2,u1,dx)==0) = 0;
-			clear dx dy	% free some memory...
-		end
-
 		fxy(isnan(fxy)) = 0;
 
-		% computes maximum absolute gradient (median-style), normalizes, saturates and duplicates in 3-D matrix
-		r = repmat(max(min(fxy/nmedian(abs(fxy),1 - scut/100),1),-1),[1,1,3]);
+		% computes maximum absolute gradient (median-style), normalizes,
+		% saturates and duplicates in 3-D matrix
+		li = 1 - abs(sind(el)); % light amplitude (experimental)
+		r = repmat(max(min(li*fxy/nmedian(abs(fxy),1 - lcut/100),1),-1),[1,1,3]);
+		rp = (1 - abs(r)).^ct;
 	
 		% applies contrast using exponent
-		rp = (1 - abs(r)).^c;
 		I = I.*rp;
 	
 		% lighter for positive gradient
-		k = find(r > 0);
-		I(k) = I(k) + (1 - rp(k));
+		I(r>0) = I(r>0) + (1 - rp(r>0));
 				
 	end
-	
-	% set novalues / NaN to black color
+
+	% set novalues / NaN to nancolor
 	[i,j] = find(isnan(z));
 	if ~isempty(i)
 		I(sub2ind(size(I),repmat(i,1,3),repmat(j,1,3),repmat(1:3,size(i,1),1))) = repmat(novalue_color,size(i,1),1);
 	end
+	
+	% lake option
+	if lake
+		klake = islake(z);
+	else
+		klake = 0;
+	end
+	
 	% set the seacolor for 0 values
 	if ~isempty(csea)
-		[i,j] = find(z==0);
+		[i,j] = find(z==0 | klake);
 		if ~isempty(i)
 			I(sub2ind(size(I),repmat(i,1,3),repmat(j,1,3),repmat(1:3,size(i,1),1))) = repmat(csea(end,:),size(i,1),1);
@@ -326,54 +643,110 @@
 	end
 
-	hh = imagesc(x,y,I);
-else
-	hh = imagesc(x,y,z);
-	colormap(cmap);
-	text(mean(x),mean(y),'SPLASH!','Color','c','FontWeight','bold','HorizontalAlignment','center')
-end
-
-orient tall
-axis xy, axis equal, axis tight
-
+	if wmark
+		I = watermark(I,wmark);
+	end
+	txt = '';
+	
+else
+	
+	I = repmat(shiftdim(sea_color,-1),size(z));
+	cmap = repmat(sea_color,[256,1]);
+	txt = 'Mak Byur!';	% Splash !
+end
+
+% -------------------------------------------------------------------------
+% --- ends the function when 'noplot' option is on
+if noplot
+	varargout{1} = struct('x',x,'y',y,'z',z,'rgb',I);
+	return
+end
+
+% -------------------------------------------------------------------------
+% --- plots the RGB image
+hh = imagesc(x,y,I);
+
+if ~isempty(txt)
+	text(mean(x),mean(y),txt,'Color',sea_color/4, ...
+		'FontWeight','bold','HorizontalAlignment','center')
+end
+
+orient tall; axis xy
+if strcmpi(axeq,'auto')
+	axis equal
+end
+axis tight
 xlim = [min(x),max(x)];
 ylim = [min(y),max(y)];
-zlim = [min(z(:)),max(z(:))];
-
-% axis basemap style
-if dec | dms
+zlim = [min([z(z(:) ~= novalue);zmin]),max([z(z(:) ~= novalue);zmax])];
+
+if dms
+	% approximates X-Y aspect ratio for this latitude (< 20-m precision for 1x1° grid)
+	xyr = cos(mean(y)*pi/180);
+else
+	xyr = 1;
+end
+
+bw0 = max(diff(xlim)*xyr,diff(ylim))/100;
+bwy = bw*bw0; % Y border width = 1%
+bwx = bwy/xyr; % border width (in degree of longitude)
+
+
+% -------------------------------------------------------------------------
+% --- Axis basemap style
+if dec || dms
 	axis off
 
-	if dms
-		% approximates X-Y aspect ratio for this latitude (< 20-m precision for 1x1° grid)
-		xyr = cos(mean(y)*pi/180);
-	else
-		xyr = 1;
-	end
-	set(gca,'DataAspectRatio',[1,xyr,1])
-	bwy = 0.008*diff(ylim); % Y border width = 1%
-	bwx = bwy/xyr; % border width (in degree of longitude)
-
-	% transparent borders
-	patch([xlim(1)-bwx,xlim(2)+bwx,xlim(2)+bwx,xlim(1)-bwx],ylim(1) - bwy*[0,0,1,1],'k','FaceColor','none','clipping','off')
-	patch([xlim(1)-bwx,xlim(2)+bwx,xlim(2)+bwx,xlim(1)-bwx],ylim(2) + bwy*[0,0,1,1],'k','FaceColor','none','clipping','off')
-	patch(xlim(1) - bwx*[0,0,1,1],[ylim(1)-bwy,ylim(2)+bwy,ylim(2)+bwy,ylim(1)-bwy],'k','FaceColor','none','clipping','off')
-	patch(xlim(2) + bwx*[0,0,1,1],[ylim(1)-bwy,ylim(2)+bwy,ylim(2)+bwy,ylim(1)-bwy],'k','FaceColor','none','clipping','off')
-
+	if strcmpi(axeq,'manual')
+		ppos = get(gcf,'PaperPosition');
+		apos = get(gca,'Position');
+		xyf = (xyr*diff(xlim)/apos(3)/ppos(3))/(diff(ylim)/apos(4)/ppos(4));
+		if xyf >= 1
+			set(gca,'Position',[apos(1),apos(2),apos(3),apos(4)/xyf]);
+		else
+			set(gca,'Position',[apos(1),apos(2),apos(3)*xyf,apos(4)]);
+		end
+	end
+	if strcmpi(axeq,'auto')
+		if diff(xlim)*xyr <= diff(ylim)
+			set(gca,'DataAspectRatio',[1,xyr,1])
+		else
+			set(gca,'DataAspectRatio',[1/xyr,1,1])
+		end
+	end
+
+	if bw > 0
+		% transparent borders
+		patch([xlim(1)-bwx,xlim(2)+bwx,xlim(2)+bwx,xlim(1)-bwx],ylim(1) - bwy*[0,0,1,1],'k','FaceColor','none','clipping','off')
+		patch([xlim(1)-bwx,xlim(2)+bwx,xlim(2)+bwx,xlim(1)-bwx],ylim(2) + bwy*[0,0,1,1],'k','FaceColor','none','clipping','off')
+		patch(xlim(1) - bwx*[0,0,1,1],[ylim(1)-bwy,ylim(2)+bwy,ylim(2)+bwy,ylim(1)-bwy],'k','FaceColor','none','clipping','off')
+		patch(xlim(2) + bwx*[0,0,1,1],[ylim(1)-bwy,ylim(2)+bwy,ylim(2)+bwy,ylim(1)-bwy],'k','FaceColor','none','clipping','off')
+	end
 	dlon = {'E','W'};
 	dlat = {'N','S'};
-
-	if dec
-		ddx = dtick(diff(xlim));
-		ddy = dtick(diff(ylim));
+	if fbold
+		fw = 'bold';
 	else
-		ddx = dtick(diff(xlim),1);
-		ddy = dtick(diff(ylim),1);
-	end
-
+		fw = 'normal';
+	end
+	
+	if ddx == 0
+		ddx = dtick(diff(xlim),dms);
+	end
+	if ddy == 0
+		ddy = dtick(diff(ylim),dms);
+	end
 	xtick = (ddx*ceil(xlim(1)/ddx)):ddx:xlim(2);
 	for xt = xtick(1:2:end)
 		dt = ddx - max(0,xt + ddx - xlim(2));
 		patch(repmat(xt + dt*[0,1,1,0]',[1,2]),[ylim(1) - bwy*[0,0,1,1];ylim(2) + bwy*[0,0,1,1]]','k','clipping','off')
-		text(xt,ylim(1) - bwy,deg2dms(xt,dlon,dec),'FontSize',fs,'HorizontalAlignment','center','VerticalAlignment','top');
+		if fs > 0
+			if ~isempty(regexp(tpos,'north','once'))
+				text(xt,ylim(2) + 1.2*bwy,deg2dms(xt,dlon,dec),'FontSize',fs,'FontWeight',fw, ...
+					'HorizontalAlignment','center','VerticalAlignment','bottom');
+			else
+				text(xt,ylim(1) - 1.2*bwy,deg2dms(xt,dlon,dec),'FontSize',fs,'FontWeight',fw, ...
+					'HorizontalAlignment','center','VerticalAlignment','top');
+			end
+		end
 	end
 
@@ -382,50 +755,77 @@
 		dt = ddy - max(0,yt + ddy - ylim(2));
 		patch([xlim(1) - bwx*[0,0,1,1];xlim(2) + bwx*[0,0,1,1]]',repmat(yt + dt*[0,1,1,0]',[1,2]),'k','clipping','off')
-		text(xlim(1) - 1.1*bwx,yt,deg2dms(yt,dlat,dec),'FontSize',fs,'HorizontalAlignment','right','VerticalAlignment','middle');
-	end
-end
-
-% scale legend
+		if fs > 0
+			if ~isempty(regexp(tpos,'east','once'))
+				text(xlim(2) + 1.2*bwx,yt,deg2dms(yt,dlat,dec),'FontSize',fs,'FontWeight',fw, ...
+					'HorizontalAlignment','center','VerticalAlignment','top','rotation',90);
+			else
+				text(xlim(1) - 1.2*bwx,yt,deg2dms(yt,dlat,dec),'FontSize',fs,'FontWeight',fw, ...
+					'HorizontalAlignment','center','VerticalAlignment','bottom','rotation',90);
+			end
+		end
+	end
+end
+
+% -------------------------------------------------------------------------
+% --- Scales legend
 if scale
 	%wsc = diff(xlim)*0.01;
-	wsc = bwx;
-	xsc = xlim(2) + wsc*4;
-
-	% elevation scale (colorbar)
-	zscale = linspace(zlim(1),zlim(2),length(cmap));
+	wsc = bw0;
+	xsc = xlim(2) + wsc*2 + bwx;
+
+	if wmark
+		cmap = watermark(cmap,wmark);
+	end
+
+	% -- elevation scale (colorbar)
+	zscale = linspace(zmin,zmax,length(cmap))';
 	yscale = linspace(0,diff(ylim)/2,length(cmap));
+	ddz = dtick(dz*max(0.5*xyr*diff(xlim)/yscale(end),1));
+	ztick = (ddz*ceil(zscale(1)/ddz)):ddz:zscale(end);
+ 	rgbscale = ind2rgb(uint16(round((zscale - zmin)*(size(cmap,1) - 1)/dz) + 1),cmap);
 	ysc = ylim(1);
-	ddz = dtick(dz*max(0.5*xyr*diff(xlim)/yscale(end),1));
-	ztick = (ddz*ceil(zlim(1)/ddz)):ddz:zlim(2);
-	patch(xsc + repmat(wsc*[-1;1;1;-1],[1,length(cmap)]), ...
-		ysc + [repmat(yscale,[2,1]);repmat(yscale + diff(yscale(1:2)),[2,1])], ...
-		repmat(zscale,[4,1]), ...
-		'EdgeColor','flat','LineWidth',.1,'FaceColor','flat','clipping','off')
-	colormap(cmap)
-	caxis([zmin,zmax])
+	hold on
+	imagesc(xsc + wsc*[-1,1]/2,ysc + yscale,repmat(rgbscale,1,2),'clipping','off');
 	patch(xsc + wsc*[-1,1,1,-1],ysc + yscale(end)*[0,0,1,1],'k','FaceColor','none','Clipping','off')
-	text(xsc + 2*wsc + zeros(size(ztick)),ysc + (ztick - zlim(1))*0.5*diff(ylim)/diff(zlim),num2str(ztick'), ...
-		'HorizontalAlignment','left','VerticalAlignment','middle','FontSize',8)
+	text(xsc + 2*wsc + zeros(size(ztick)),ysc + (ztick - zscale(1))*0.5*diff(ylim)/diff(zscale([1,end])),num2str(ztick'), ...
+		'HorizontalAlignment','left','VerticalAlignment','middle','FontSize',6)
+	% indicates min and max Z values
+	text(xsc,ysc - bwy/2,sprintf('%g m',roundsd(zlim(1),3)),'FontWeight','bold', ...
+		'HorizontalAlignment','left','VerticalAlignment','top','FontSize',6)
+	text(xsc,ysc + .5*diff(ylim) + bwy/2,sprintf('%g m',roundsd(zlim(2),3)),'FontWeight','bold', ...
+		'HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',6)
 	
-	% distance scale (in case of DMS only)
+	% frees axes only if not hold on
+	if ~holdon
+		hold off
+	end
+	
+	% -- distance scale (in km)
 	if dms
-		degkm = 6370*pi/180;
-		dkm = dtick(diff(ylim)*degkm);
-		ysc = ylim(2) - 0.5*dkm/degkm;
-		patch(xsc + wsc*[-1,-1,0,0],ysc + dkm*0.5*[-1,1,1,-1]/degkm,'k','FaceColor',grey,'clipping','off')
-		if dkm > 1
-			skm = sprintf('%g km',dkm);
-		else
-			skm = sprintf('%g m',dkm*1000);
-		end
-		text(xsc,ysc,skm,'rotation',-90,'HorizontalAlignment','center','VerticalAlignment','bottom', ...
-			'Color',grey,'FontWeight','bold')
-	end
-
+		fsc = degkm;
+	else
+		fsc = zratio/1e3;
+	end
+	dkm = dtick(diff(ylim)*fsc);
+	ysc = ylim(2) - 0.5*dkm/fsc;
+	patch(xsc + wsc*[-1,-1,0,0],ysc + dkm*0.5*[-1,1,1,-1]/fsc,'k','FaceColor',grey,'clipping','off')
+	if dkm > 1
+		skm = sprintf('%g km',dkm);
+	else
+		skm = sprintf('%g m',dkm*1000);
+	end
+	text(xsc,ysc,skm,'rotation',-90,'HorizontalAlignment','center','VerticalAlignment','bottom', ...
+			'Color',grey,'FontWeight','bold','FontSize',6)
 end
 
 
 if nargout > 0
-	h = hh;
+	varargout{1} = hh;
+end
+if nargout > 1
+	varargout{2} = I;
+end
+if nargout > 2
+	varargout{3} = z;
 end
 
@@ -455,8 +855,8 @@
 end
 
-if deg & dlim <= 2/60
+if deg && dlim <= 2/60
 	% less than 2 minutes: base 36
 	m = 10^floor(log10(dlim*36))/36;
-elseif deg & dlim <= 2
+elseif deg && dlim <= 2
 	% less than 2 degrees: base 6
 	m = 10^floor(log10(dlim*6))/6;
@@ -485,5 +885,6 @@
 else
 	xa = abs(x) + 1/360000;
-	sd = sprintf('%d%c',floor(xa),176);	% ASCII char 176 is the degree sign
+	%sd = sprintf('%d%c',floor(xa),176);	% ASCII char 176 is the degree sign
+	sd = sprintf('%d°',floor(xa));
 	sm = '';
 	ss = '';
@@ -502,9 +903,10 @@
 end
 
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function z = fillgap(x,y,z)
 % GRIDDATA is not efficient for large arrays, but has great advantage to be
-% included in Matlab core functions! To optimize interpolation, we
-% reduce the number of relevant data by building a mask of surrounding
+% included in Matlab's core functions! To optimize interpolation, we
+% reduce the amount of relevant data by building a mask of all surrounding
 % pixels of novalue areas... playing with linear index!
 
@@ -529,2 +931,151 @@
 [i,j] = ind2sub(sz,k);
 k2 = sub2ind(fliplr(sz),j,i); % switched i and j: k2 is linear index in row order
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function k = islake(z)
+% ISLAKE mask of zero gradient on 3x3 tiles
+% We use diff matrix in row and column directions, and shift it to build
+% a single vectorized test of surrounding pixels. To do this we must
+% concatenate unit vectors in different combinations...
+
+dx = diff(z,1,2);	% differences in X direction
+dy = diff(z,1,1);	% differences in Y direction
+u1 = ones(size(z,1),1);	% row unit vector 
+u2 = ones(1,size(z,2));	% column unit vector
+u2r = u2(2:end);
+
+% index of the tiles center pixel
+k = ( ...
+	[u2;dy] == 0 & [dy;u2] == 0 & ...
+	[u1,dx] == 0 & [dx,u1] == 0 & ...
+	[u1,[dx(2:end,:);u2r]] == 0 & [[dx(2:end,:);u2r],u1] == 0 & ...
+	[u1,[u2r;dx(1:end-1,:)]] == 0 & [[u2r;dx(1:end-1,:)],u1] == 0 ...
+);
+
+% now extends it to surrounding pixels
+k(1:end-1,:) = (k(1:end-1,:) | k(2:end,:));
+k(2:end,:) = (k(2:end,:) | k(1:end-1,:));
+k(:,1:end-1) = (k(:,1:end-1) | k(:,2:end));
+k(:,2:end) = (k(:,2:end) | k(:,1:end-1));
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function s = isrgb(x,n)
+
+if nargin < 2
+	n = 0;
+end
+if isnumeric(x) && (n == 1 && all(size(x) == [1,3]) || n == 0 && size(x,2) == 3) ...
+		&& all(x(:) >= 0 & x(:) <= 1)
+	s = 1;
+else
+	s = 0;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function s = isperc(x)
+
+if isnumeric(x) && isscalar(x) && x >= 0 && x <= 100
+	s = 1;
+else
+	s = 0;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function s = isvec(x,n)
+
+if nargin < 2
+	n = 2;
+end
+if isnumeric(x) && numel(x) == n
+	s = 1;
+else
+	s = 0;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function y=roundsd(x,n)
+
+og = 10.^(floor(log10(abs(x)) - n + 1));
+y = round(x./og).*og;
+y(x==0) = 0;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function y = watermark(x,n)
+
+if nargin < 2
+	n = 2;
+end
+
+if n == 0
+    y = x;
+else
+    y = (x/n + 1 - 1/n);
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [s,v] = checkparam(arg,nam,func,val)
+
+switch func2str(func)
+	case 'isscalar'
+		num = 1;
+		mes = 'scalar value';
+	case 'isperc'
+		num = 1;
+		mes = 'percentage scalar value';
+	case 'isvec'
+		num = 1;
+		if nargin < 4
+			val = 2;
+		end
+		mes = sprintf('%d-element vector',val);
+	case 'isrgb'
+		num = 1;
+		mes = '[R,G,B] vector with 0.0 to 1.0 values';
+	case 'ischar'
+		num = 0;
+		mes = 'string';
+		if nargin > 3
+			mes = sprintf('%s (%s)',mes,strjoin(val,' or '));
+		end
+	otherwise
+		num = 1;
+		mes = 'value';
+end
+
+s = 0;
+v = [];
+k = find(strcmpi(arg,nam));
+if ~isempty(k)
+	if (k + 1) <= length(arg) ...
+			&& (~num || isnumeric(arg{k+1})) ...
+			&& (nargin < 4 && func(arg{k+1}) ...
+				|| (nargin > 3 && (strcmp(func2str(func),'ischar') && ismember(arg{k+1},val)) ...
+					 || strcmp(func2str(func),'isvec') && func(arg{k+1},val)))
+		v = arg{k+1};
+		s = 1;
+	else
+		error('%s option must be followed by a valid %s.',upper(nam),mes)
+	end
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function s=strjoin(c,d)
+%STRJOIN Join cell array of strings
+%(this is for Matlab versions < 2013a backward compatibility)
+
+if nargin < 2
+	d = '';
+end
+n = numel(c);
+ss = cell(2,n);
+ss(1,:) = reshape(c,1,n);
+ss(2,1:n-1) = {d};
+s = [ss{:}];
Index: /issm/trunk-jpl/externalpackages/dem/license.txt
===================================================================
--- /issm/trunk-jpl/externalpackages/dem/license.txt	(revision 23662)
+++ /issm/trunk-jpl/externalpackages/dem/license.txt	(revision 23663)
@@ -1,24 +1,24 @@
-Copyright (c) 2013, François Beauducel
+Copyright (c) 2016, François Beauducel
 All rights reserved.
 
-Redistribution and use in source and binary forms, with or without 
-modification, are permitted provided that the following conditions are 
+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 
+    * 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 
+    * 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 
+
+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.