use matlab to draw a curving column

已知轴线方程：

思想：把圆柱按轴线分若干段，在每个断面画圆，圆的法线就轴线的切向方向。柱面画若干平行于轴线的线条，其中用到平行单位矢量点积等于1，不平行的结果是[-1，1）。

程序如下：

% rotation
clc;close all;clear;
t=30;
p=60;
t=t/180*pi;
p=p/180*pi;
T1=[1,     0,     0;
    0,cos(p),-sin(p);
    0,sin(p),cos(p)];
X=[0;0;1];  % normal
T2=[sin(t),-cos(t),0;
    cos(t), sin(t),0;
    0,           0,1];
Rx=T1*T2*X; % retated normal
N=100;      % the number of sample point
t=linspace(0,pi*2,N); % sample points
P=[sin(t);cos(t);zeros(1,N)]; % points on circle
Rp=T1*T2*P; % retated circle
plot3(Rp(1,:),Rp(2,:),Rp(3,:),'r');hold on;
plot3(P(1,:),P(2,:),P(3,:),'b');
axis([-1,1,-1,1,-1,1]);
plot3([0,X(1)],[0,X(2)],[0,X(3)],'b','linewidth',2);
plot3([0,Rx(1)],[0,Rx(2)],[0,Rx(3)],'r','linewidth',2);
view([-67,14])
figure;
t=linspace(0,pi*2,N);
Q=[cos(t*2);sin(t*2);t];  % the axes equation of column
plot3(Q(1,:),Q(2,:),Q(3,:),'r');
T=T1*T2;
Nq=[Q(:,2:end)-Q(:,1:end-1)];Nq=[Nq,Nq(:,end)];
R=0.2; % radii
P=R*P;
hold on;
M=7;
m=round(linspace(0,N,M+1));m(1)=[];

for k=1:length(Q);
    q=Q(:,k);
    q=q/sqrt(sum(q.^2));
    p=acos(q(3));
    t=angle(q(1)+i*q(2));
    T1=[1,     0,     0;
        0,cos(p),-sin(p);
        0,sin(p),cos(p)];
    T2=[sin(t),-cos(t),0;
        cos(t), sin(t),0;
        0,           0,1];
    T=T1*T2;  % T is rotation matrix in 3D
    Rp=T*P;
    plot3(Q(1,k)+Rp(1,:),Q(2,k)+Rp(2,:),Q(3,k)+Rp(3,:),'k');
    if k==1;
        Bx(1:M,k)=Q(1,k)+[Rp(1,m(1:M))]';
        By(1:M,k)=Q(2,k)+[Rp(2,m(1:M))]';
        Bz(1:M,k)=Q(3,k)+[Rp(3,m(1:M))]';
    else
        nn=Q(:,k)-Q(:,k-1);nn=nn/sqrt(sum(nn.^2));
        Nn=repmat(nn,1,N);
        for n=1:M;          
            B1=Q(1,k)+Rp(1,:)-Bx(n,k-1);
            B2=Q(2,k)+Rp(2,:)-By(n,k-1);
            B3=Q(3,k)+Rp(3,:)-Bz(n,k-1);
            NN=[B1;B2;B3];
            NNs=repmat(sqrt(sum(NN.^2)),3,1);
            NN=NN./NNs;
            S=sum(Nn.*NN);
            [s,po]=max(S)
            Bx(n,k)=Q(1,k)+Rp(1,po);
            By(n,k)=Q(2,k)+Rp(2,po);
            Bz(n,k)=Q(3,k)+Rp(3,po);
        end
    end
end
for k=1:M;
    plot3(Bx(k,:),By(k,:),Bz(k,:),'k');
end
axis equal;

结果如下：