function [x1,y1]=cubicRefine(x0,y0) % Wenjie He, July 05
% function [x1,y1]=cubicRefine(x0,y0) Do one refinement based on cubic splines. 
% Input is knot positions and values
% Output is new positions and values for one level of refinement

s1=max(size(x0));
s2=max(size(y0));
if s1 ~= s2
   error('The sizes of input do not match.');
end
if min(size(x0)) > 1
   error('The first argument must be a vector.');
end
if min(size(y0)) > 1
   error('The second argument must be a vector.');
end
if s1 < 4
   error('The size of the input is too small, s>=4.');
end

x1=zeros(1,2*s1-3);
y1=zeros(1,2*s1-3);

x1(1)=x0(1);
y1(1)=y0(1);
x1(2*s1-3)=x0(s1);
y1(2*s1-3)=y0(s1);
x1(2)=(x0(1)+x0(2))/2;
y1(2)=(y0(1)+y0(2))/2;
x1(2*s1-4)=(x0(s1-1)+x0(s1))/2;
y1(2*s1-4)=(y0(s1-1)+y0(s1))/2;

if s1 == 4
   x1(3)=(x0(2)+x0(3))/2;
   y1(3)=(y0(2)+y0(3))/2;
else
   x1(3)=(3*x0(2)+x0(3))/4;% knots near boundary are special.
   y1(3)=(3*y0(2)+y0(3))/4;
   x1(2*s1-5)=(x0(s1-2)+3*x0(s1-1))/4;
   y1(2*s1-5)=(y0(s1-2)+3*y0(s1-1))/4;
end

if s1 == 5
   x1(4)=(3*x0(2)+10*x0(3)+3*x0(4))/16;
   y1(4)=(3*y0(2)+10*y0(3)+3*y0(4))/16;
end

if s1 >= 6
   x1(4)=(3*x0(2)+11*x0(3)+2*x0(4))/16;
   y1(4)=(3*y0(2)+11*y0(3)+2*y0(4))/16;
   for i=3:s1-3
      x1(2*i-1)=(x0(i)+x0(i+1))/2;
      y1(2*i-1)=(y0(i)+y0(i+1))/2;
	   x1(2*i)=(x0(i)+6*x0(i+1)+x0(i+2))/8;
      y1(2*i)=(y0(i)+6*y0(i+1)+y0(i+2))/8;
   end
   x1(2*s1-6)=(2*x0(s1-3)+11*x0(s1-2)+3*x0(s1-1))/16;
   y1(2*s1-6)=(2*y0(s1-3)+11*y0(s1-2)+3*y0(s1-1))/16;
end