%
% Script File: h5curvesEa.m = curves1224.m
% Bezier Curve Fit on the example at bottom of P.317 modified
%                     from Mathew and Fink
% (1) Bezier Curve 
% (2) Cubic B-Spline
% (3) Cubic Spline Interpolant
%     By using N (N=7 here) control points
%
% X0=[0, 1, 2, 3, 4, 5, 6];
% Y0=[3, 5, 1, 3, 5, 2, 4];
fin=fopen('dataEa.txt','r'); fgetl(fin);
N=7; % N= size(X0) =size(Y0);
X0=fscanf(fin,'%f',[1,N]);
Y0=fscanf(fin,'%f',[1,N]);
n=41;
U=linspace(0,1,n);
V=1.0-U;
%
% (1) Bezier Curve
%
Xz=[];  Yz=[];
for i=1:(N-1)/3
  j=3*(i-1)+1;
  Xj=V.^3*X0(j)+3*(V.^2).*U*X0(j+1)+3*V.*(U.^2)*X0(j+2)+U.^3*X0(j+3);
  Yj=V.^3*Y0(j)+3*(V.^2).*U*Y0(j+1)+3*V.*(U.^2)*Y0(j+2)+U.^3*Y0(j+3);
  Xz=[Xz,Xj];
  Yz=[Yz,Yj];
end
%
% (2) Cubic B-Spline
%
Bx=[]; By=[];
for k=2:N-2
  Bk=(1-U).^3/6*X0(k-1)+(U.^3/2-U.^2+2/3)*X0(k)+...
       (-U.^3+U.^2+U+1/3)/2*X0(k+1)+U.^3/6*X0(k+2);
  Ck=(1-U).^3/6*Y0(k-1)+(U.^3/2-U.^2+2/3)*Y0(k)+...
     (-U.^3+U.^2+U+1/3)/2*Y0(k+1)+U.^3/6*Y0(k+2);
  Bx=[Bx,Bk];
  By=[By,Ck];
end
%
% (3) Cubic Spline Interpolant
%
Xs=0:0.1:6;
Ys=spline(X0,Y0,Xs);
plot(X0,Y0,'ro',Xz,Yz,'c-',Bx,By,'m-',Xs,Ys,'b-'); axis([0,6, 0,6])
legend('7 Points','Bezier','B-Spline','Cubic Spline','Location','Southeast');
title('Cubic Bezier Curve, B-Spline, Cubic Spline Interpolant from 7 Control Points')