function windschief

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

deg=pi/180;
close('all')
P1=[4; 2; 1];
v1=[-1; 6; 5];
P2=[-10; -4; 3];
v2=[-2; -1; 8];
format compact
format long
n=cross(v1,v2)
Nenner=n'*n
t1Stern=det([P2-P1,v2,n])/Nenner
t2Stern=det([P2-P1,v1,n])/Nenner
P1Stern=P1+t1Stern*v1
P2Stern=P2+t2Stern*v2
t1Gra=[min([0 t1Stern])-3/norm(v1) max([0 t1Stern])+3/norm(v1)];
t2Gra=[min([0 t2Stern])-3/norm(v2) max([0 t2Stern])+3/norm(v2)];
plot3(P1(1)+t1Gra*v1(1),P1(2)+t1Gra*v1(2),P1(3)+t1Gra*v1(3),'k',...
      P2(1)+t2Gra*v2(1),P2(2)+t2Gra*v2(2),P2(3)+t2Gra*v2(3),'k',...
      [P1Stern(1) P2Stern(1)],[P1Stern(2) P2Stern(2)],[P1Stern(3) P2Stern(3)],'--k',...
      [P1(1) P2(1)],[P1(2) P2(2)],[P1(3) P2(3)],'ok','MarkerSize',10,'linewidth',1.5);
hold('on')
plot3([P1Stern(1) P2Stern(1)],[P1Stern(2) P2Stern(2)],[P1Stern(3) P2Stern(3)],'.k','MarkerSize',30)
d=P2Stern-P1Stern;
d=d/norm(d);
M1=(P1+P1Stern)/2;
M2=(P2+P2Stern)/2;
A=M1+0.5*d-1.75*v1/norm(v1);
B=M1+0.5*d+1.75*v1/norm(v1);
plot3([A(1) B(1)],[A(2) B(2)],[A(3) B(3)],'k')
Spitze(B,v1,1.0,30*deg)
A=M2+0.5*d-1.75*v2/norm(v2);
B=M2+0.5*d+1.75*v2/norm(v2);
plot3([A(1) B(1)],[A(2) B(2)],[A(3) B(3)],'k')
Spitze(B,v2,1.0,30*deg)
axis('equal');
axis([-9 5 -5 0 -8 4])
grid('on');
set(gca,'xticklabel',[]);
set(gca,'yticklabel',[]);
set(gca,'zticklabel',[]);
Schrift(P1,' {\itP}_0','VerticalAlignment','baseline')
Schrift(P2,' {\itQ}_0')
Schrift(P1Stern,' {\itP}({\itt}*)','VerticalAlignment','baseline')
Schrift(P2Stern,' {\itQ}({\its}*)')
Schrift(M1,' {\itP}({\itt})','VerticalAlignment','baseline')
Schrift(M2,' {\itQ}({\its})')
Schrift(M1+d,'{\it\bfv}','HorizontalAlignment','right','VerticalAlignment','middle')
Schrift(M2+0.75*d,'{\it\bfw}','HorizontalAlignment','right','VerticalAlignment','middle')
print('-deps','windschief')

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

function Spitze(P,v,L,phi);
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');
  end
end
