function Kegel

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

global ax
ax=[inf -inf inf -inf inf -inf];
clc;
close('all')
c=inline('3+sin(3*t).*cos(5*t)./(1+sin(exp(t)).^2)','t');
k=inline('3','t');
t=linspace(-pi,pi,48);
[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:2);6];

hold('on')
Keg(M,S,c,length(t))

M(1)=M(1)+7;
S(1)=S(1)+14;
Keg(M,S,c,length(t))

M(1)=M(1)+8;
S(1)=M(1);
Keg(M,S,k,length(t))

M(1)=M(1)+6.7;
S(1)=S(1)+14;
Keg(M,S,k,length(t))

t=linspace(-pi,pi,6);
[x,y]=pol2cart(t,c(t));
x=x+M(1)+8.3;
xx=x([2:end 1]);
yy=y([2:end 1]);
M0=[[(x.^2+x.*xx+xx.^2)*(yy-y)';(y.^2+y.*yy+yy.^2)*(x-xx)']/(6*polyarea(x,y));0];
M=M0+[0;M(2)-M0(2);0];
S=[M(1:2);S(3)];
Keg(M,S,c,length(t))

M(1)=M(1)+6.5;
S(1)=M(1)+7;
Keg(M,S,c,length(t))

hold('off')
axis('equal')
view(0,30)
axis([ax(1:5) 4.5])
axis('off')
print('-deps','Kegel','-zbuffer')

function Keg(M,Sp,c,n)
global ax
v=Sp-M;
t=linspace(0,2*pi,n);
s=[0 1];
[T,S]=meshgrid(t,s);
X=M(1)+S*v(1)+(1-S).*c(T).*cos(T);
Y=M(2)+S*v(2)+(1-S).*c(T).*sin(T);
Z=M(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])
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
