function f = sBA ( x,y,z, m, n )
% sBA - BA algorithm, fits a B-spline surface through the given points
%	f = sBA ( x,y,z, xmax, ymax ) or f = sBA ( [x y z], [xmax, ymax] )
%	x,y,z - vectors with fitting points, it is assumed that z=f(x,y)
%	xmax, ymax - size of the output spline coeffitient matrix
%	f - matrix with fitted surface B-spline control points
%
% see also: sMBA

L = 1; U = 1; % Left/Upper offsets from the first element in the matrix

delta = zeros(m+U+3,n+L+2);
w2    = zeros(m+U+3,n+L+2);
% in i direction, matrices are extended by U from upper, by 3 from lower side
% in j direction, matrices are extended by L from left , by 2 from right side

for c = 1:length(x),
if( round(c/100) == c/100 ) disp (sprintf('    %d %%', round(100*c/length(x)))); end;
	i = floor(x(c))-1;  j = floor(y(c))-1;
        I = i+(0:3)+1+U; J = j+(0:3)+1+L;
	s = x(c)-floor(x(c));  t = y(c)-floor(y(c));
        S = 1/6 * [(1-s)^3  (3*s^3-6*s^2+4) (-3*s^3+3*s^2+3*s+1)  s^3]' ;
        T = 1/6 * [(1-t)^3  (3*t^3-6*t^2+4) (-3*t^3+3*t^2+3*t+1)  t^3]' ;
	w = S*T'; W2=sum(sum(w.^2));

	fi= w*z(c)/W2;
	delta(I,J) = delta(I,J) + w.^2.*fi;
	w2   (I,J) = w2   (I,J) + w.^2;
end;

% revolving 0 and 360 degrees meridians 
delta((1:2)+U, : ) = delta((1:2)+U, : ) + delta((m+1:m+2)+U, : );
w2   ((1:2)+U, : ) = w2   ((1:2)+U, : ) + w2   ((m+1:m+2)+U, : );
delta(m+U,:) = delta(m+U,:) + delta(1+U,:);
w2   (m+U,:) = w2   (m+U,:) + w2   (1+U,:);

% overlapping parallels from the opposite sides of the globe accross poles
if m == 1,   t = 1; else, t = [ (m/2+1:m) (1:m/2) ]; end
delta((1:m)+U,(1:2)+L)   = delta((1:m)+U, (1:2)   +L) ...
		         + delta(  t  +U, (1:-1:0)+L);
delta((1:m)+U,(n-1:n)+L) = delta((1:m)+U, (n-1:n) +L) ...
		         + delta(  t  +U, (n+1:-1:n)+L);
w2   ((1:m)+U,(1:2)+L)   = w2   ((1:m)+U, (1:2)   +L) ...
		         + w2   (  t  +U, (1:-1:0)+L);
w2   ((1:m)+U,(n-1:n)+L) = w2   ((1:m)+U, (n-1:n) +L) ...
		         + w2   (  t  +U, (n+1:-1:n)+L);

f = (w2~=0).*delta ./ ( w2 + (w2==0) );  % f = delta./w2, for w2~=0
% returning central portion of the matrix
f = f((1:m)+U,(1:n)+L);

return;