function x = nice(v,fac);
% NICE  'nice' scaling vector.
%
%  X = NICE(V), where V is a 2-element vector [A,B], returns
%      an arbitrarily-but-nicely spaced vector X such that [A,B] 
%      is surrounded by the limits of X.
%
%  X = NICE(V,F) , where F is an integer, increases the resolution
%      of X by a factor F.
%
%      V may also be a NxM matrix; in this case, only its minimum 
%      and maximum value are used.
%      If minimum and maximum of V are equal, a warning message is
%      given.
%
%  Example:
%	           x = nice( [0.72 2.39] )
%	           y = nice( peaks )
%  returns
%	           x == [ 0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0  2.2  2.4 ] ,
%	           y == [ -8  -6  -4  -2  0  2  4  6  8  10 ] .
%
%  (c) 08-02-1995  Wolfgang Erasmi, Dept. Marine Physics, IfM Kiel, Germany
%      werasmi@ifm.uni-kiel.d400.de
%
%  Note: this version uses max9 and min9 for data that may contain NaN
%        values; replace them by max, min if necessary.
%	Revised  31-08-1995  W. Erasmi
%           Changes to make function work even for max == min 
%           and size(V) == [1 1].

% constants
mlbd1 = 4;		% boundary for 1st increase of no.s of levels
mfac1 = 2;		%     increasing factor for this case
mlbd2 = 1.8;	% boundary for 2nd increase of no.s of levels
mfac2 = 5;		%     increasing factor for this case

% error handling & defaults
if nargin < 2,
	mult = 1;		% default multiplier
else
	mult = max(fix(fac),1);
end

a = min9(min9(v));    			% minimum
b = max9(max9(v));  			% maximum

if ( a == b ),
%	fprintf(1,'%c',7);	% beep
	disp(' ');
	disp('nice: Warning: Upper equals lower boundary.');
	if a == 0,
		a = -1;	b = 1;
	else
		a = a - abs(a);
		b = b + abs(b);
	end
end

d = b-a;							% difference
od = round(log10(d)-0.5);	% order of magnitude of d; this looks somewhat
									% weird but has to be used instead of fix(log10(d))
									% because of numerical inaccuracy in log10.
scal = 10^od;
ds = d/scal;					% relative delta (between 0.1 and 1)

% calculate stepwidth of vector; use stepwise increase of resolution,
%                                if necessary
if ds > mlbd1,
	dlev = 1/mult * scal;
elseif ds > mlbd2,
	dlev = 1/(mult*mfac1) * scal;
else 
	dlev = 1/(mult*mfac2) * scal;
end

% create output vector
lo = floor(a/dlev)*dlev;	
hi =  ceil(b/dlev)*dlev;
x = (lo:dlev:hi);		
