function CavalieriZylinder

% http://geometrie.hs-nb.de/MATLAB-Programme/CavalieriZylinder.m   Bearbeitungsstand: 6.10.2004
% (C) Martin Nitschke, nitschke@hs-nb.de

global ax
ax=[inf -inf inf -inf inf -inf];
clc;
close('all')
c=inline('3+2*sin(3*t).*cos(5*t)./(1+sin(exp(t)).^2)','t');
nSch=12;
t=linspace(-pi,pi,32);
[x,y]=pol2cart(t,c(t));
xx=x([2:end 1]);
yy=y([2:end 1]);
M=[[(x.^2+x.*xx+xx.^2)*(yy-y)';(y.^2+y.*yy+yy.^2)*(x-xx)']/(6*polyarea(x,y));0];
S=[M(1)+10;M(2);6];
v=S-M;
v=v/norm(v);
hold('on')
Zyl(M,S,c,length(t))
M=M+[12;0;0];
S=S+[12;0;0];
h=linspace(0,norm(S-M),nSch);
for k=2:length(h)
  Zyl(M+h(k-1)*v,M+h(k-1)*v+[0;0;(S(3)-M(3))/nSch],c,length(t))
end
M=M+[17;0;0];
S=S+[17;0;0];
S(1)=M(1);
v=S-M;
v=v/norm(v);
h=linspace(0,norm(S-M),nSch);
for k=2:length(h)
  Zyl(M+h(k-1)*v,M+h(k)*v,c,length(t))
end
M=M+[8;0;0];
S=S+[8;0;0];
Mx=M(1);
Zyl(M,S,c,length(t))

c=inline('repmat(3,size(t))','t');
nSch=12;
t=linspace(-pi,pi,24);
[x,y]=pol2cart(t,c(t));
xx=x([2:end 1]);
yy=y([2:end 1]);
M=[[(x.^2+x.*xx+xx.^2)*(yy-y)';(y.^2+y.*yy+yy.^2)*(x-xx)']/(6*polyarea(x,y));0];
S=[M(1)+10;M(2);6];
v=S-M;
v=v/norm(v);
M=M+[Mx+12;0;0];
S=S+[Mx+12;0;0];
Zyl(M,S,c,length(t))
M=M+[11;0;0];
S=S+[11;0;0];
h=linspace(0,norm(S-M),nSch);
for k=2:length(h)
  Zyl(M+h(k-1)*v,M+h(k-1)*v+[0;0;(S(3)-M(3))/nSch],c,length(t))
end
M=M+[16;0;0];
S=S+[16;0;0];
S(1)=M(1);
v=S-M;
v=v/norm(v);
h=linspace(0,norm(S-M),nSch);
for k=2:length(h)
  Zyl(M+h(k-1)*v,M+h(k)*v,c,length(t))
end
M=M+[7;0;0];
S=S+[7;0;0];
Zyl(M,S,c,length(t))
hold('off')
axis('equal')
axis('off')
view(0,20)
axis(ax)
print('-deps','CavalieriZylinder','-opengl')

function Zyl(M1,M2,c,n)
global ax
v=M2-M1;
t=linspace(0,2*pi,n);
s=[0 1];
[T,S]=meshgrid(t,s);
X=M1(1)+S*v(1)+c(T).*cos(T);
Y=M1(2)+S*v(2)+c(T).*sin(T);
Z=M1(3)+S*v(3);
surf(X,Y,Z,T)
c=0.5+0.375*(sin(pi*[0:length(t)-1]/(length(t)-1))).^2;
colormap(repmat(c',1,3));
fill3(X(1,:),Y(1,:),Z(1,:),mean(c)*[1 1 1])
fill3(X(end,:),Y(end,:),Z(end,:),mean(c)*[1 1 1])
if min(X(:))<ax(1)
  ax(1)=min(X(:));
end
if max(X(:))>ax(2)
  ax(2)=max(X(:));
end
if min(Y(:))<ax(3)
  ax(3)=min(Y(:));
end
if max(Y(:))>ax(4)
  ax(4)=max(Y(:));
end
if min(Z(:))<ax(5)
  ax(5)=min(Z(:));
end
if max(Z(:))>ax(6)
  ax(6)=max(Z(:));
end
