function  Psi = sMBA( x, y, z, maxn, epsilon, makeplot )
% sMBA - MBA algorithm, fitting surface with splines (algorithm by Prof. Wolberg)
%	   M = sMBA( x, y, z, maxn, epsilon ) 
%      or  M = sMBA( [x y z], maxn, epsilon )
%	x,y,z - vectors of scatterd points. It is assumed that z = f(x,y) 
%	maxn - maximal size of the mesh, (DEFAULT 16)
%	epsilon - bound for maximal error (DEFAULT 1/100*mean(z) )
%       makeplot - string indicating if surface should be plotted
%	   'p' means plotting, anything else means no plotting
%
%       M - n*(n+1) matrix with spline node coefficients M(fi,teta)
%	   M(1,1) corresponds to fi=-pi, teta=-pi/2
%          fi-space is in 2*pi/n wide segments (-pi  :2*pi/n: pi  ) 
%          teta-space is in pi/n wide segments (-pi/2 : pi/n: pi/2)
%	       convenient for use with for use with sPlotSurfM
%
% see also: sBA, sFit, sPlotSurfM
 
if ( min(size(x))==3 ),
   if nargin >= 4, makeplot = maxn; else makeplot='n'; end;
   if nargin >= 2, maxn = y; else maxn=64; end;
   y=x(:,2); z=x(:,3); x=x(:,1);
   if nargin >= 3, epsilon = z; else epsilon= 1/100 * mean( z ); end;
else
   if nargin < 6, makeplot='n'; end;
   if nargin < 5, epsilon= 1/100 * mean( z ); end;
   if nargin < 4, maxn=64; end;
end;

if makeplot(1) == 'p' | makeplot(1) == 'P', makeplot = 'p'; end; 
maxerror = epsilon+1;


n0 = 2;
Psi = mean(z) * ones(n0, n0+1);	   % initial spline control points
z = z - mean(z);		   % subtract initial spline approximation

n = size(Psi,1);
%disp('relative max error |  mean error');     % banner for displaying error through iterations

% maxn = 64;
if makeplot == 'p',
   f = gcf;  % get the current figure number
end;

while  n <= maxn/2, % & maxerror > epsilon
disp(n);
       n = n*2;			       %double the resolution
       xn = n*x; yn = n*y;	       %scaling (0:1)->(0:n)   
       Fi = sBA( xn, yn, z, n,n+1);    %approximation of z
       z  = z - sValue(Fi,xn,yn);      %decrease z by the approximation 
       Psi= sRefine(Psi) + Fi;	       %double the number of control points and
				       % increase them by the approximation
       maxerror = max(abs(z));   
       % The rest is just plotting    
       if makeplot == 'p',
	   disp([maxerror/epsilon mean(abs(z))]);	   % diplaying error
	   % f = gcf;
	   if(n>2), 
	       figure(f); whitebg('w'); sSurfM( Psi); axis('off');
	       % xlabel('el'); ylabel('az'); title(sprintf('Mesh size %dx%d', n,n));
	       % grid on; 
	   end; 
	   figure(f+1); whitebg('w'); sPlotSurf((Psi)); axis('off');
	   % grid on; 
	   view(3);drawnow; 
	   figure(f);
	   f = f+2;
       end; %plotting
end;

return