function Schraubkoerper

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

R=6; h=12;
deg=pi/180;
syms('t','s','u','v','real');
clc;
close('all')
c1=0.7*(3+sin(3*t)*cos(5*t)/(1+sin(exp(t))^2))*cos(t);
c2=0.9*(3+sin(3*t)*cos(5*t)/(1+sin(exp(t))^2))*sin(t);
ti=[0:5:360]*deg;
nSchr=2.5;
Dm=[0 nSchr]*2*pi;
m=[R*[cos(s);sin(s)];h*s/(2*pi)];
mInl=inline(['[' char(m(1)) ',' char(m(2)) ',' char(m(3)) '];'],'s');
si=linspace(Dm(1),Dm(2),nSchr*length(ti))';
mi=mInl(si);
Phi=VRK(c1,c2,t,m,s,u,v);
xRKInl=inline(vectorize(Phi(1)),'u','v');
yRKInl=inline(vectorize(Phi(2)),'u','v');
zRKInl=inline(vectorize(Phi(3)),'u','v');
[U,V]=meshgrid(ti,si);
X=xRKInl(U,V);
Y=yRKInl(U,V);
Z=zRKInl(U,V);
surf(X,Y,Z,U);
c=0.5+0.375*(sin(pi*[0:size(V,2)-1]/(size(V,2)-1))).^2;
colormap(repmat(c',1,3));
hold('on')
plot3(X([1 end],:)',Y([1 end],:)',Z([1 end],:)','k','linewidth',2.0)
Br=unique(round(linspace(1,size(X,2),32)));
plot3(X(:,Br),Y(:,Br),Z(:,Br),'k')
PM=Polygon3M([mi(:,1) mi(:,2) mi(:,3)]',1.0,'ok','erase','none','MarkerSize',5);
zRg=[min(Z(:));max(Z(:))]+(max(Z(:))-min(Z(:)))*[-0.2;0.1];
plot3([0;0],[0;0],zRg,'k','linewidth',2.0)
plot3(mi(:,1),mi(:,2),repmat(zRg(1)/2,size(si)),'k')
sRef=1;
Spitze([mi(sRef,1:2)';zRg(1)/2], [mi(sRef,1:2)';0],1.5,30*deg);
Spitze([0;0;zRg(1)/2]          ,-[mi(sRef,1:2)';0],1.5,30*deg);
plot3([mi(sRef,1) 0],[mi(sRef,2) 0],zRg([1 1])/2,'k')
Schrift([mi(sRef,1:2)'/2;zRg(1)/2],'\itR','HorizontalAlignment','center','VerticalAlignment','top','erase','none')
[mx mxi]=min(mi(1:round(end/nSchr),1));
mRef=mi(mxi,:)';
Spitze(mRef,        [0;0;-1],2.5,30*deg,'Erase','none');
Spitze(mRef+[0;0;h],[0;0; 1],2.5,30*deg,'Erase','none');
plot3(mRef([1 1]),mRef([2 2]),mRef(3)+[0 h],'k','Erase','none')
Schrift(mRef+[0;0;h/2],'\ith_{ }','HorizontalAlignment','right','VerticalAlignment','middle','erase','none')

[U,V]=meshgrid(ti,[0 h*nSchr]);
o=15+1.5*R;
surf(o+R*cos(U),R*sin(U),V,U)
drawPolygon3(@plot3,[o+1.05*PM(1,:);1.05*PM(2,:);PM(3,:)],'ok','MarkerSize',5);
plot3(o+[0;0],[0;0],zRg,'k','linewidth',2.0)
plot3(o+mi(:,1),mi(:,2),repmat(zRg(1)/2,size(si)),'k')
plot3(o+[mi(sRef,1) 0],[mi(sRef,2) 0],zRg([1 1])/2,'k')
Spitze([[o;0]+mi(sRef,1:2)';zRg(1)/2], [mi(sRef,1:2)';0],1.5,30*deg,'erase','none');
Spitze([o;0;zRg(1)/2]          ,-[mi(sRef,1:2)';0],1.5,30*deg,'erase','none');
Schrift([[o;0]+mi(sRef,1:2)'/2;zRg(1)/2],'\itR','HorizontalAlignment','center','VerticalAlignment','top','erase','none')

o=o+1.5*R;
surf(o+R*U,zeros(size(V))-R,V,U)
drawPolygon3(@plot3,[o+mod(2*pi*R*PM(3,:)/h,2*pi*R);zeros(size(PM(3,:)))-1.05*R;PM(3,:)],'ok','MarkerSize',5);
for k=1:floor(nSchr)
  drawPolygon3(@plot3,[o+[0 2*pi*R];[0 0]-R;h*[k k]],'k','erase','none');
  Schrift([o;-R;h*(k-0.5)],'\ith','HorizontalAlignment','right','VerticalAlignment','middle','erase','none')
end
Schrift([o+pi*R;-R;0],'{\itU}=2\pi{\itR}','HorizontalAlignment','center','VerticalAlignment','top','erase','none')
view(-11,10)
axis('equal')
axis([min(X(:)) max(X(:))+o+2*pi*R min(Y(:)) max(Y(:)) zRg(1) 0.87*max(Z(:))])
axis('off')
shading('interp')
print('-deps','Schraubkoerper','-opengl')

function Phi=VRK(c1,c2,t,m,s,u,v);
dm=simple(diff(m,s));
T=simple(dm/sqrt(dot(dm,dm)));
Z=sym([0;0;1]);
N1=cross(Z,T);
N1=simple(N1/sqrt(dot(N1,N1)));
N2=cross(N1,T);
N2=simple(N2/sqrt(dot(N2,N2)));
Phi=subs(m,s,v)+subs(c1,t,u)*subs(N1,s,v)+subs(c2,t,u)*subs(N2,s,v);

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

function drawPolygon3(draw,P,varargin);
feval(draw,P(1,:),P(2,:),P(3,:),varargin{:})

function Q=Polygon3M(P,d,varargin)
lastK=1;
Q=P(:,1);
drawPolygon3(@plot3,P(:,1),varargin{:});
for k=2:size(P,2)
  v=P(:,k)-P(:,k-1);
  t=linspace(0,1,round(norm(v)/d));
  if length(t)==1
    pl=norm(P(:,lastK)-P(:,k))>d;
  else
    pl=1;
  end
  if pl
    lastK=k;
    R=repmat(P(:,k-1),1,length(t))+repmat(t,length(P(:,k-1)),1).*repmat(v,1,length(t));
    drawPolygon3(@plot3,R,varargin{:})
    Q=[Q R];
  end
end

function Spitze(P,v,L,phi,varargin);
deg=pi/180;
if nargin<=2
  L=1;
end
if nargin<=3
  phi=7.5*deg;
end
v=L*v./repmat(sqrt(dot(v,v)),size(P,1),1);
Q=P-v;
t=[0:6:360]*deg;
for k=1:size(P,2)
  N=tan(phi/2)*null(v(:,k)');
  for w=2:length(t)
    plot3(...
      [P(1,k) Q(1,k)+cos(t(w))*L(:,k)*N(1,1)+sin(t(w))*L(:,k)*N(1,2) Q(1,k)+cos(t(w-1))*L(:,k)*N(1,1)+sin(t(w-1))*L(:,k)*N(1,2) P(1,k)],...
      [P(2,k) Q(2,k)+cos(t(w))*L(:,k)*N(2,1)+sin(t(w))*L(:,k)*N(2,2) Q(2,k)+cos(t(w-1))*L(:,k)*N(2,1)+sin(t(w-1))*L(:,k)*N(2,2) P(2,k)],...
      [P(3,k) Q(3,k)+cos(t(w))*L(:,k)*N(3,1)+sin(t(w))*L(:,k)*N(3,2) Q(3,k)+cos(t(w-1))*L(:,k)*N(3,1)+sin(t(w-1))*L(:,k)*N(3,2) P(3,k)],...
      'k',varargin{:});
    end
  end
