function Konoid

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

clc;
close('all')
deg=pi/180;
global ax
syms('t','u','v','real');
n=[0;1;0];                     % Normalenvektor der Richtebene
P=[0;0;2]; w=[0;1;0];          % Leitgerade
c=[(3+sin(5*t)/7)*cos(t);(2+sin(7*t)/5)*sin(t);0];       % Leitkurve
Phi=Kon(P,w,c,t,n,u,v);
xKonInl=inline(vectorize(Phi(1)),'u','v');
yKonInl=inline(vectorize(Phi(2)),'u','v');
zKonInl=inline(vectorize(Phi(3)),'u','v');
ti=[0:6:360]*deg;
[U,V]=meshgrid(ti,[0 1]);
X=xKonInl(U,V);
Y=yKonInl(U,V);
Z=zKonInl(U,V);
surf(X,Y,Z,U);
c=0.5+0.375*(sin(pi*[0:size(U,2)-1]/(size(U,2)-1))).^2;
colormap(repmat(c',1,3));
xRgKon=[min(X(:)) max(X(:))];
yRgKon=[min(Y(:)) max(Y(:))];
zRgKon=[min(Z(:)) max(Z(:))];
ax=[xRgKon yRgKon yRgKon];
A=[mean(xRgKon);2*yRgKon(1)-1*yRgKon(2);zRgKon(1)];
B=double(SchnittGerEb(P,w,A,n));
a=cross(B-A,n);
a=0.5*(xRgKon(2)-xRgKon(1))*a/norm(a);
hold('on')
drawPolygon(@plot3,[X(1,:);Y(1,:);Z(1,:)],'k','linewidth',2.0)
vRef=round(length(ti)/2);
Schrift([X(1,vRef);Y(1,vRef);Z(1,vRef)],{'Leit-' 'kurve'},'HorizontalAlignment','right','VerticalAlignment','top','erase','none')
S=[[X(end,:);Y(end,:);Z(end,:)] B];
tS=(S-repmat(P,1,size(S,2)))'*w;
tRg=[min(tS) max(tS)]+(max(tS)-min(tS))*[-0.4 0.1];
drawPolygon(@plot3,[P+tRg(1)*w P+tRg(2)*w],'k')
Schrift(P+tRg(2)*w,{'Leit-' 'gerade'},'HorizontalAlignment','right','VerticalAlignment','middle','erase','none')
B=1.5*B-0.5*A;
drawPolygon(@fill3,[A-a A+a B+a B-a],mean(c)*[1 1 1],'EdgeColor','none')
Schrift(B-a,{'Richt-' 'ebene'},'HorizontalAlignment','right','VerticalAlignment','bottom','erase','none')
hold('off')
view(-40,36)
axis('equal')
axis([ax(1) ax(2) 0.8*ax(3) ax(4) 0*ax(5) 0.5*ax(6)])
axis('off')
print('-deps','Konoid','-opengl');

function Phi=Kon(P,w,c,t,n,u,v)
k=subs(c,t,u);  
Phi=simple(k+v*(SchnittGerEb(P,w,k,n)-k));

function S=SchnittGerEb(A,v,B,n)
syms('t','real');
S=simple(A+solve((A+t*v-B)'*n,t)*v);

function drawPolygon(draw,P,varargin);
global ax
feval(draw,P(1,:),P(2,:),P(3,:),varargin{:})
if min(P(1,:))<ax(1)
  ax(1)=min(P(1,:));
end
if max(P(1,:))>ax(2)
  ax(2)=max(P(1,:));
end
if min(P(2,:))<ax(3)
  ax(3)=min(P(2,:));
end
if max(P(2,:))>ax(4)
  ax(4)=max(P(2,:));
end
if min(P(3,:))<ax(5)
  ax(5)=min(P(3,:));
end
if max(P(3,:))>ax(6)
  ax(6)=max(P(3,:));
end

function Schrift(P,txt,varargin)
text(P(1),P(2),P(3),txt,'FontName','TimesNewRoman','FontSize',24,varargin{:});
