%clear
%数据点
x=1:8;
y=0.5*x.^2+3;
% y=x;
qq0=[sqrt(2);sqrt(2)];qqe=[sqrt(64/65);sqrt(1/65)];
% qqe=qq0;
q=[x;y];
q=[
34.4507 64.8979 92.7423
53.2133 62.9103 78.7443
72.6174 58.5389 64.6679
90.4922 53.8726 51.7908
107.6262 50.2527 40.0300
124.3115 48.5382 29.7099
142.9174 47.8914 20.1672
165.0593 46.1233 11.5884
187.7633 41.4142 5.4394
210.3851 32.3890 1.3837
232.5009 19.6438 -0.7049
253.2947 4.2544 -1.1379
]';
qq0=[-0.6489 0.7079 0.2789]';%给定首数据点切矢
qqe=[0.1227 0.7740 -0.6212]';%给定末数据点切矢
[~,m]=size(q);
m=m-1;
n=m+2;
%求节点矢量
L=zeros(1,m);
for i=1:m
L(i)=sqrt((q(1,i+1)-q(1,i))^2+(q(2,i+1)-q(2,i))^2+(q(3,i+1)-q(3,i))^2);%%%维度
end
Ls=sum(L);
U=zeros(1,n+5);
U(1,n:end)=1;
for i=1:m
U(3+i+1)=sum(L(1:i))/Ls;
end
%求控制点
delta=U(2:end)-U(1:end-1);
delta=[delta,0];
a=zeros(1,n-1);b=a;c=a;e=[a;a;a]; %%%维度
% for i=3:n-1
% a(i-1)=delta(i+2)^2/(delta(i)+delta(i+1)+delta(i+2));
% b(i-1)=delta(i+2)*(delta(i)+delta(i+1))/(delta(i)+delta(i+1)+delta(i+2))+delta(i+1)*(delta(i+2)+delta(i+3))/(delta(i+1)+delta(i+2)+delta(i+3));
% c(i-1)=delta(i+1)^2/(delta(i+1)+delta(i+2)+delta(i+3));
% e(:,i-1)=(delta(i+1)+delta(i+2)).*q(:,i);
% end
for i=2:n-2%a,b,c,e的值
a(i)=delta(i+3)^2/(delta(i+1)+delta(i+2)+delta(i+3));
b(i)=delta(i+3)*(delta(i+1)+delta(i+2))/(delta(i+1)+delta(i+2)+...
delta(i+3))+delta(i+2)*(delta(i+3)+delta(i+4))/(delta(i+2)+...
delta(i+3)+delta(i+4));
c(i)=delta(i+2)^2/(delta(i+2)+delta(i+3)+delta(i+4));
e(:,i)=(delta(i+2)+delta(i+3)).*q(:,i);
end
% %边界条件——切矢条件
% a(1)=1;
% b(1)=0;
% c(1)=0;
% e(:,1)=q(:,1)+qq0/3*delta(4);
% a(n-1)=0;
% b(n-1)=0;
% c(n-1)=1;
% e(:,n-1)=q(:,m+1)-(delta(n+1)/3)*qqe;
%边界条件——自由端点条件
a(1)=2-delta(4)*delta(5)/(delta(4)+delta(5))^2;
b(1)=delta(4)/(delta(4)+delta(5))*(delta(5)/(delta(4)+delta(5))-delta(4)/(delta(4)+delta(5)+delta(6)));
c(1)=delta(4)^2/(delta(4)+delta(5))/(delta(4)+delta(5)+delta(6));
e(:,1)=q(:,1)+q(:,2);
a(n-1)=-delta(n+1)^2/(delta(n)+delta(n+1))*(delta(n-1)+delta(n)+delta(n+1));
b(n-1)=delta(n+1)/(delta(n)+delta(n+1))*(delta(n+1)/(delta(n-1)+delta(n)+delta(n+1))-delta(n)/(delta(n)+delta(n+1)));
c(n-1)=delta(n)*delta(n+1)/(delta(n)+delta(n+1))^2-2;
e(:,n-1)=-q(:,m+1)-q(:,m); %%m
% % %边界条件——抛物线条件
% a(1)=1-delta(4)*delta(5)/(delta(4)+delta(5))^2;
% b(1)=delta(4)/(delta(4)+delta(5))*(delta(5)/(delta(4)+delta(5))-delta(4)/(delta(4)+delta(5)+delta(6)));
% c(1)=delta(4)^2/(delta(4)+delta(5))/(delta(4)+delta(5)+delta(6));
% e(:,1)=q(:,1)/3+2*q(:,2)/3;
% a(n-1)=-delta(n+1)^2/(delta(n)+delta(n+1))*(delta(n-1)+delta(n)+delta(n+1));
% b(n-1)=delta(n+1)/(delta(n)+delta(n+1))*(delta(n+1)/(delta(n-1)+delta(n)+delta(n+1))-delta(n)/(delta(n)+delta(n+1)));
% c(n-1)=delta(n)*delta(n+1)/(delta(n)+delta(n+1))^2-2;
% e(:,n-1)=-q(:,m+1)/3-2*q(:,m)/3; %%m
%矩阵方程
myA=zeros(n-1,n-1);
myA(1,1:3)=[a(1),b(1),c(1)];
myA(n-1,n-3:n-1)=[a(n-1),b(n-1),c(n-1)];
for i=2:n-2
myA(i,i-1:i+1)=[a(i),b(i),c(i)];
end
d=(myA\e')';
d=[q(:,1),d,q(:,end)];
u=0:0.01:1;
p=[u;u;u]; %%%维度
for i=1:length(u)
p(:,i)=[0;0;0]; %%%维度
for j=1:n+1
p(:,i)=p(:,i)+Bbasis(j-1,3,u(i),U)*d(:,j);
end
end
p(:,end)=q(:,end);
% plot(p(1,1:end),p(2,1:end),'r')
% hold on
% plot(d(1,:),d(2,:),'b')
% hold on
% plot(q(1,:),q(2,:),'.')
plot3(p(1,:),p(2,:),p(3,:),'k')
hold on
plot3(d(1,:),d(2,:),d(3,:),'b-o')
hold on
plot3(q(1,:),q(2,:),q(3,:),'*r')
