function Kegelschnitte

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

global ax AbbNr
ax=[inf -inf inf -inf inf -inf]; 
AbbNr=0;
[SEll,FEll]=KS('[-1.1;-1;-1]+s*[0;2;0]+t*[2;0;1.3]','k','linewidth',5.0);
[SPar,FPar]=KS('[-0.6;-1;-1]+s*[0;2;0]+t*[2;0;2]','k','linewidth',5.0);
[SHyp,FHyp]=KS('[-0.6;-1;-1]+s*[0;2;0]+t*[1.5;0;2]','k','linewidth',5.0);
KS('[-1;-1;-0.25]+s*[0;2;0]+t*[2;0;0.5]','.k','MarkerSize',60,'erase','none');
KS('[-1;-1;-1]+s*[0;2;0]+t*[2;0;2]','k','linewidth',5.0);
KS('[-0.75;-1;-1]+s*[0;2;0]+t*[1.5;0;2]','k','linewidth',5.0);
Ell(SEll,FEll);
Par(SPar,FPar);
Hyp(SHyp,FHyp);

function [S,F]=KS(Eb,varargin)
global ax AbbNr
close all
AbbNr=AbbNr+1;
deg=pi/180;
gKegel=inline('x^2+y^2-z^2','x','y','z');
PEbene=inline(Eb,'s','t');
name=['Kegelschnitte-' int2str(AbbNr)];
figure('numbertitle','off','name',name);
hold('on')
drawPolygon(@plot3,[PEbene(0,0) PEbene(1,0) PEbene(1,1) PEbene(0,1) PEbene(0,0)],'k','linewidth',2.5);
syms('s','t','real')
v1=double(diff(eval(Eb),s));
v2=double(diff(eval(Eb),t));
P0=double(eval(Eb)-s*v1-t*v2);
f=solve((P0+s*v1+t*v2)'*[1;0;0],(P0+s*v1+t*v2)'*[0;1;0]);
F=P0+double(f.s)*v1+double(f.t)*v2;
P=PEbene(s,t);
sSchn=solve(simple(gKegel(P(1),P(2),P(3))),s);
if sSchn(1)==sSchn(2)
  tS=[];
else
  tS=double(solve(sSchn(1)-sSchn(2),t));
  if length(tS)==1
    S=PEbene(subs(sSchn(1),t,tS),tS);
  else
    S1=PEbene(subs(sSchn(1),t,tS(1)),tS(1));
    S2=PEbene(subs(sSchn(2),t,tS(2)),tS(2));
    if norm(S1-F)<norm(S2-F)
      F1=F;
      F2=S2+S1-F1;
    else
      F2=F;
      F1=S1+S2-F2;
    end
    S=[S1 S2];
    F=[F1 F2];
  end
end
for k=1:length(sSchn)
  sInl=inline(vectorize(sSchn(k)),'t');
  [tP,sP]=fplot(sInl,[0 1],5e-4);
  tP=unique([tP;tS]);
  sP=sInl(tP);
  epsIm=find(abs(imag(sP))<=sqrt(eps));
  sP(epsIm)=real(sP(epsIm));
  sP(find(imag(sP)))=NaN;
  if length(sP)==1
    sP=repmat(sP,size(tP));
  end
  P=zeros(3,length(tP));
  for p=1:length(tP)
    P(:,p)=PEbene(sP(p),tP(p));
  end
  drawPolygon(@plot3,P,varargin{:})
end
h=linspace(-1,1,32);
lambda=[0:360]*deg;
[H,Lambda]=meshgrid(h,lambda);
surf(H.*cos(Lambda),H.*sin(Lambda),H,Lambda)
c=0.5+0.375*(sin(pi*[0:size(H,1)-1]/(size(H,1)-1))).^2;
colormap(repmat(c',1,3));
shading('interp')
axis(ax)
axis('equal')
axis('off')
view(-66,18);
axis([-1 1 -1 1 -1 1])
print('-deps','-loose',name,'-opengl');
view(0,0);
Schrift([1;0;-0.97],'\alpha  ','FontSize',96,'erase','none','HorizontalAlignment','right','VerticalAlignment','baseline')
t=[0:45]*deg;
r=0.7;
drawPolygon(@plot3,[1-r*cos(t);zeros(size(t));-1+r*sin(t)],'k','erase','none');
axis([-1 1 -1 1 -1 1])
hold('off')
print('-deps','-loose',[name 'o'],'-opengl');

function Ell(S,F)
global AbbNr
close all
h=0.1;
a=norm(S(:,2)-S(:,1))/2;
e=a-norm(S(:,1)-F(:,1));
b=sqrt(a^2-e^2);
ax=[-a-1.5*h a+2*h -b-0.5*h b+1.5*h];
AbbNr=AbbNr+1;
deg=pi/180;
name=['Kegelschnitte-' int2str(AbbNr)];

figure('numbertitle','off','name',name);
[xG,yG]=meshgrid([ax(1) -a -e e a ax(2)],[ax(3) -b b ax(4)]);
t=[0:360]*deg;
plot(xG,yG,xG',yG','color',0.5*[1 1 1]);
hold('on')
plot(a*cos(t),b*sin(t),'k','linewidth',2.0)
plot(ax(1:2),[0 0],'k')
plot([0 0],ax(3:4),'k')
Schrift([-a;0],'{\itS}_2','HorizontalAlignment','right','VerticalAlignment','top')
Schrift([-e;0],'{\itF}_2','HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ e;0],'{\itF}_1','HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ a;0],'{\itS}_1','HorizontalAlignment','left','VerticalAlignment','top')
Schrift([0; b],'{\itN}_1','HorizontalAlignment','right','VerticalAlignment','top')
Schrift([0;-b],'{\itN}_2','HorizontalAlignment','right','VerticalAlignment','bottom')
plot([-e 0 e],[0 b 0],'k')
Schrift([e;b]/2,'\ita','rotation',atan(b/e)/deg-90,'HorizontalAlignment','center','VerticalAlignment','bottom')
xRef=-(1.5*e-0.5*a);
yRef=b*sqrt(1-(xRef/a)^2);
Schrift([xRef;yRef],'{\itP}','color',0.5*[1 1 1],'HorizontalAlignment','center','VerticalAlignment','baseline')
plot([-e xRef e],[0 yRef 0],'color',0.5*[1 1 1])
Schrift([xRef+e;yRef]/2,'_{ }\itr_1','color',0.5*[1 1 1],'FontSize',30,'rotation',atan(yRef/(xRef-e))/deg,...
  'HorizontalAlignment','left','VerticalAlignment','top')
Schrift([xRef-e;yRef]/2,'\itr_2','color',0.5*[1 1 1],'FontSize',30,'rotation',atan(yRef/(xRef+e))/deg,...
  'HorizontalAlignment','left','VerticalAlignment','bottom')
Spitze([ax(2);0],[1;0],h,30*deg,'k')
Schrift([ax(2);0.2*h],'\itx','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Spitze([0;ax(4)],[0;1],h,30*deg,'k')
Schrift([0;ax(4)+0.4*h],'\ity ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Schrift([-a;ax(3)],'\it-a','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([-e;ax(3)],'\it-e','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ 0;ax(3)],'0','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ e;ax(3)],'\ite','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ a;ax(3)],'\ita','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ax(1);-b],'\it-b ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','middle')
Schrift([ax(1); 0],'0 ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','middle')
Schrift([ax(1); b],'\itb ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','middle')
axis('equal')
axis(ax)
axis('off')
hold('off')
print('-depsc',name);

figure('numbertitle','off','name',[name 'R']);
plot(a*cos(t),b*sin(t),'k','linewidth',2.0)
hold('on')
phi=[0:15:345]*deg;
p=b^2/a;
ex=e/a;
r=p./(1+ex*cos(phi));
[x,y]=pol2cart(phi,r);
for k=1:length(x)
  plot([e e+x(k) -e],[0 y(k) 0],'k')
end
axis('equal')
axis(ax)
axis('off')
hold('off')
print('-deps',[name 'R']);

function Par(S,F)
global AbbNr
close all
h=0.1;
p=2*norm(S-F);
ax=[-3*p 6*p 0 0];
ax(3:4)=[-1 1].*sqrt(2*p*ax(2));
AbbNr=AbbNr+1;
deg=pi/180;
name=['Kegelschnitte-' int2str(AbbNr)];

figure('numbertitle','off','name',name);
[xG,yG]=meshgrid([ax(1) -p/2 p/2 ax(2)],[ax(3) ax(4)]);
t=linspace(ax(3),ax(4),100);
plot(xG,yG,xG',yG','color',0.5*[1 1 1]);
hold('on')
plot(t.^2/(2*p),t,'k','linewidth',2.0)
plot(ax(1:2),[0 0],'k')
plot([0 0],ax(3:4),'k')
Schrift([0;0],'{\itS}','HorizontalAlignment','right','VerticalAlignment','top')
Schrift([p/2;0],'{\itF}','HorizontalAlignment','center','VerticalAlignment','top')
xRef=p;
yRef=sqrt(2*p*xRef);
plot([-p/2 -p/2],ax(3:4),'color',0.5*[1 1 1],'linewidth',1.0);
plot([-p/2 xRef p/2],[yRef yRef 0],'color',0.5*[1 1 1]);
Schrift([xRef;yRef],'{\itP}','color',0.5*[1 1 1],'HorizontalAlignment','center','VerticalAlignment','baseline')
Schrift([-p/2;yRef/2],'{\itl}_{ }','color',0.5*[1 1 1],'HorizontalAlignment','right','VerticalAlignment','middle')
Schrift([xRef/2;yRef],'{\itr}_{ }','color',0.5*[1 1 1],'HorizontalAlignment','right','VerticalAlignment','baseline')
Schrift([(xRef+p/2)/2;yRef/2],'{\itr}','color',0.5*[1 1 1],'HorizontalAlignment','left','VerticalAlignment','middle')

Spitze([ax(2);0],[1;0],h,30*deg,'k')
Schrift([ax(2);0.2*h],'\itx','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Spitze([0;ax(4)],[0;1],h,30*deg,'k')
Schrift([0;ax(4)+0.4*h],'\ity ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Schrift([-p/2;ax(3)],'\it-p/2','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Schrift([ 0;ax(3)],'0','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ p/2;ax(3)],'\itp/2','FontSize',30,'HorizontalAlignment','left','VerticalAlignment','top')
Schrift([ax(1); 0],'0 ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','middle')
axis('equal')
axis(ax)
axis('off')
hold('off')
print('-depsc',name);

figure('numbertitle','off','name',[name 'R']);
plot(t.^2/(2*p),t,'k','linewidth',2.0)
hold('on')
phi=[0:15:345]*deg;
r=p./(1+cos(phi));
[x,y]=pol2cart(phi,r);
for k=1:length(x)
  plot([p/2 p/2-x(k) ax(2)],[0 y(k) y(k)],'k')
end
axis('equal')
axis(ax)
axis('off')
hold('off')
print('-deps','-loose',[name 'R']);

function Hyp(S,F)
global AbbNr
close all
h=0.1;
a=norm(S(:,2)-S(:,1))/2;
e=a+norm(S(:,1)-F(:,1));
b=sqrt(e^2-a^2);
ax=[-1.6*e 1.6*e 0 0];
ax(3:4)=[-1 1].*b*sqrt((ax(1)/a)^2-1);
AbbNr=AbbNr+1;
deg=pi/180;

name=['Kegelschnitte-' int2str(AbbNr)];
figure('numbertitle','off','name',name);
[xG,yG]=meshgrid([ax(1) -e -a a e ax(2)],[ax(3) -b b ax(4)]);
t=linspace(-acosh(ax(2)/a),acosh(ax(2)/a),100);
plot(xG,yG,xG',yG','color',0.5*[1 1 1]);
hold('on')
plot( a*cosh(t), b*sinh(t),'k','linewidth',2.0)
plot(-a*cosh(t), b*sinh(t),'k','linewidth',2.0)
plot(ax(1:2),[0 0],'k')
plot([0 0],ax(3:4),'k')
Schrift([-a;0],'{\itS}_2','HorizontalAlignment','left','VerticalAlignment','top')
Schrift([-e;0],'{\itF}_2','HorizontalAlignment','right','VerticalAlignment','top')
Schrift([ e;0],'{\itF}_1','HorizontalAlignment','left','VerticalAlignment','top')
Schrift([ a;0],'{\itS}_1','HorizontalAlignment','right','VerticalAlignment','top')
m=b/a;
plot(ax(3:4)/m, ax(3:4),'k')
plot(ax(3:4)/m,-ax(3:4),'k')
xRef=-(1.5*e-0.5*a);
yRef=b*sqrt((xRef/a)^2-1);
Schrift([xRef;yRef],'{\itP}','color',0.5*[1 1 1],'HorizontalAlignment','left','VerticalAlignment','baseline')
plot([-e xRef e],[0 yRef 0],'color',0.5*[1 1 1])
Schrift([xRef+e;yRef]/2,'\itr_1','color',0.5*[1 1 1],'FontSize',30,'rotation',atan(yRef/(xRef-e))/deg,...
  'HorizontalAlignment','right','VerticalAlignment','bottom')
Schrift([xRef-e;yRef]/2,'\itr_2','color',0.5*[1 1 1],'FontSize',30,'rotation',atan(yRef/(xRef+e))/deg,...
  'HorizontalAlignment','center','VerticalAlignment','top')
Spitze([ax(2);0],[1;0],h,30*deg,'k')
Schrift([ax(2);0.2*h],'\itx','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Spitze([0;ax(4)],[0;1],h,30*deg,'k')
Schrift([0;ax(4)+0.4*h],'\ity ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Schrift([-a;ax(3)],'\it-a','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([-e;ax(3)],'\it-e','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','top')
Schrift([ 0;ax(3)],'0','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ e;ax(3)],'\ite','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ a;ax(3)],'\ita','FontSize',30,'HorizontalAlignment','center','VerticalAlignment','top')
Schrift([ax(1);-b],'\it-b ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','middle')
Schrift([ax(1); 0],'0 ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','middle')
Schrift([ax(1); b],'\itb ','FontSize',30,'HorizontalAlignment','right','VerticalAlignment','middle')
axis('equal')
axis(ax)
axis('off')
hold('off')
print('-depsc',name);

figure('numbertitle','off','name',[name 'R']);
plot( a*cosh(t),b*sinh(t),'k','linewidth',2.0)
hold('on')
plot(-a*cosh(t),b*sinh(t),'k','linewidth',2.0)
phi=[0:15:345]*deg;
p=b^2/a;
ex=e/a;
r=p./(1+ex*cos(phi));
[x,y]=pol2cart(phi,r);
for k=1:length(x)
  if r(k)>0
    plot([e e-x(k) -e],[0 y(k) 0],'k')
    t=asinh(y(k)/b);
    N=[a*sinh(t);b*cosh(t)];
    N=[-N(2);N(1)]/norm(N);
    q=0.05;
    plot(e-x(k)+[-q q]*N(1),y(k)+[-q q]*N(2),'k')
  end
end
axis('equal')
axis(ax)
axis('off')
hold('off')
print('-deps',[name 'R']);

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)
if length(P)==2
  P=[P;0];
end
text(P(1),P(2),P(3),txt,'FontName','TimesNewRoman','FontSize',36,varargin{:});

function Spitze(P,v,L,phi,varargin);
if nargin<=2
  L=1;
end
if nargin<=3
  phi=7.5*pi/180;
end
ta=tan(phi/2);
v=L*v./repmat(sqrt(dot(v,v)),size(P,1),1);
n=[-v(2,:);v(1,:)];
Q=P-v+ta*n;
R=P-v-ta*n;
for k=1:size(P,2);
  fill([P(1,k) Q(1,k) R(1,k)],[P(2,k) Q(2,k) R(2,k)],varargin{:})
end
