function y = bspline3(x, n, h)

% bspline3 - cubic b-spline
% ----------------------------------------------
% y = bspline3(x, n, h)
%
% Input:
% ------
%  x - input parameter vector
%  n - offset (number of indexes)
%  h - width parameter
%
% Output:
% -------
%  y - spline function values

%--------------------------------
% Author: Matt Robbins
%--------------------------------
% $Revision$
% $Date$
%--------------------------------

%--
% handle input
%--

if nargin < 3 || isempty(h)
	h = 1;
end

if nargin < 2 || isempty(n)
	n = 0;
end

x = x - h*n;

%--
% split x into pieces
%--

x_1 = x(x < -2*h);

x_2 = x(x >= -2*h & x < -h) + 2*h;

x_3 = x(x >= -h & x < 0) + h;

x_4 = h - x(x >= 0 & x < h);

x_5 = 2*h -x(x >= h & x < 2*h);

x_6 = x(x >= 2*h);

%--
% compute y by segment
%--

y_1 = zeros(size(x_1)); y_6 = zeros(size(x_6));

y_2 = x_2.^3;

y_3 = h^3 + 3*h^2*x_3 + 3*h*x_3.^2 - 3*x_3.^3; 

y_4 = h^3 + 3*h^2*x_4 + 3*h*x_4.^2 - 3*x_4.^3; 

y_5 = x_5.^3;

%--
% assemble y
%--

y = [y_1, y_2, y_3, y_4, y_5, y_6] / h^3;
