function F2O

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

global AbbNr
AbbNr=0;
close('all');
syms a b c lambda phi u v real;
Fl([0;u;v],'u=linspace(-10,10,2);','v=linspace(-9,9,2);',[-10 10 -6 6 -9 9],'Eb');
Fl([-a a;u u;v v],'u=linspace(-10,10,2);','v=linspace(-9,9,2);',[-10 10 -6 6 -9 9],'Eb');
Fl([u u;b*u/a -b*u/a;v v],'u=linspace(-10,10,3);','v=linspace(-9,9,3);',[-10 10 -6 6 -9 9]);
Fl([a*cos(lambda);b*sin(lambda);v],'lambda=[0:pi/180:2*pi];','v=linspace(-8,8,25);',[-9 9 -4 4 -6 6]);
Fl([u^2/(2*a);u;v],'u=9*[-1:0.01:1];','v=linspace(-4,4,25);',[-9 9 -4 4 -6 6]);
Fl([a*cosh(u) -a*cosh(u);b*sinh(u) b*sinh(u);v v],'u=2*[-1:0.01:1];','v=linspace(-6,6,25);',[-9 9 -4 4 -6 6]);
Fl([a*v*cos(lambda);b*v*sin(lambda);c*v],'lambda=[0:pi/180:2*pi];','v=linspace(-6,6,25);',[-25 30 -18 18 -20 20]);
Fl([a*cos(lambda)*cos(phi);b*sin(lambda)*cos(phi);c*sin(phi)],'lambda=[0:pi/180:2*pi];','phi=[-pi/2:pi/180:pi/2];',[-4 5 -3 3 -1 1]);
Fl([a*v*cos(lambda);b*v*sin(lambda);v^2],'lambda=[0:pi/180:2*pi];','v=[0:0.02:2];',[-9 9 -6 4  2 3]);
Fl([u;v;((u/a)^2-(v/b)^2)/2],'u=linspace(-20,20,25);','v=linspace(-12,12,25);',[-20 20 -12  12 -7  1],'none');
Fl([u;v;((u/a)^2-(v/b)^2)/2],'u=linspace(-20,20,25);','v=linspace(-12,12,25);',[-20 20 -12  12 -7  1],'HP1');
Fl([u;v;((u/a)^2-(v/b)^2)/2],'u=linspace(-20,20,25);','v=linspace(-12,12,25);',[-20 20 -12  12 -7  1],'HP2');
Fl([a*cos(lambda)*cosh(phi);b*sin(lambda)*cosh(phi);c*sinh(phi)],'lambda=[0:pi/180:2*pi];','phi=[-1:0.01:1];',[-7 7 -4 4 -5 5]);
Fl([a*cos(lambda)*cosh(phi);b*sin(lambda)*cosh(phi);c*sinh(phi)],'lambda=[0:pi/180:2*pi];','phi=[-1:0.01:1];',[-7 7 -4 4 -5 5],'EH');
Fl([a*cos(lambda)*sinh(phi)  a*cos(lambda)*sinh(phi)
    b*sin(lambda)*sinh(phi)  b*sin(lambda)*sinh(phi)
    c*cosh(phi)             -c*cosh(phi)],'lambda=[0:pi/180:2*pi];','phi=[0:0.01:2];',[-15 18 -10 9 -15 13]);
  
function Fl(P,p1,p2,ax,opt)
global AbbNr
AbbNr=AbbNr+1;
if nargin<5
  opt='';
end
name=['F2O-' int2str(AbbNr)];
figure('numbertitle','off','name',name);
hold('on')
for k=1:size(P,2)
  [X,Y,Z]=FlTeil(P(:,k),p1,p2,opt);
end
view(157,44);
hold('off')
axis('equal')
if ~isempty(ax)
  axis(ax); 
end
axis('off')
print('-depsc',name,'-opengl','-loose')

function [X,Y,Z]=FlTeil(P,p1,p2,opt)
syms a b c real;
abc={5 3 4};
P=subs(P,{a b c},abc);
v1=char(symvar(p1));
v2=char(symvar(p2));
xInl=inline(vectorize(char(P(1))),v1,v2);
yInl=inline(vectorize(char(P(2))),v1,v2);
zInl=inline(vectorize(char(P(3))),v1,v2);
eval(p1);
eval(p2);
eval(['sRg=' v1 ';']);
eval(['tRg=' v2 ';']);
[S,T]=meshgrid(sRg,tRg);
X=xInl(S,T);
Y=yInl(S,T);
Z=zInl(S,T);
if max(size(X))==1
  X=repmat(X,size(T));
end
if max(size(Y))==1
  Y=repmat(Y,size(T));
end
if max(size(Z))==1
  Z=repmat(Z,size(T));
end
surf(X,Y,Z,S);
cmap=0.375+0.375*(sin(pi*[0:size(S,1)-1]/(size(S,1)-1))).^2;
colormap(repmat(cmap',1,3));
shading('interp')

switch opt
  case 'HP1'  
    a=abc{1};
    b=abc{2};
    plot3([min(X(:)) max(X(:))], (b/a)*[min(X(:)) max(X(:))],[0 0],'color',0.9*[1 1 1],'linewidth',4.0)
    plot3([min(X(:)) max(X(:))],-(b/a)*[min(X(:)) max(X(:))],[0 0],'color',0.9*[1 1 1],'linewidth',4.0)
    for h=linspace(max(Z(:))/4,max(Z(:)),5)
      tmin=acosh(max(X(:))/(sqrt(2*h)*a));
      t=linspace(-tmin,tmin,100);
      plot3( a*cosh(t)*sqrt(2*h), b*sinh(t)*sqrt(2*h),repmat(h,size(t)),'color',0.9*[1 1 1],'linewidth',4.0);
      plot3(-a*cosh(t)*sqrt(2*h), b*sinh(t)*sqrt(2*h),repmat(h,size(t)),'color',0.9*[1 1 1],'linewidth',4.0);
      plot3( a*sinh(t)*sqrt(2*h), b*cosh(t)*sqrt(2*h),repmat(-h,size(t)),'color',0.9*[1 1 1],'linewidth',4.0);
      plot3( a*sinh(t)*sqrt(2*h),-b*cosh(t)*sqrt(2*h),repmat(-h,size(t)),'color',0.9*[1 1 1],'linewidth',4.0);
    end
  case 'HP2'
    P=subs('[a*(u+v)/2;b*(u-v)/2;u*v/2]',{a b c},abc);
    xInl=inline(vectorize(char(P(1))),'u','v');
    yInl=inline(vectorize(char(P(2))),'u','v');
    zInl=inline(vectorize(char(P(3))),'u','v');
    a=abc{1};
    b=abc{2};
    plot3([min(X(:)) max(X(:))], (b/a)*[min(X(:)) max(X(:))],[0 0],'color',0.9*[1 1 1],'linewidth',4.0)
    plot3([min(X(:)) max(X(:))],-(b/a)*[min(X(:)) max(X(:))],[0 0],'color',0.9*[1 1 1],'linewidth',4.0)
    for u=linspace(min(S(:))/a+min(T(:))/b,max(S(:))/a+max(T(:))/b,13)
      plotInBox([xInl(u,[0 1]);yInl(u,[0 1]);zInl(u,[0 1])],min(S(:)),max(S(:)),min(T(:)),max(T(:)));
    end
    for v=linspace(min(S(:))/a-max(T(:))/b,max(S(:))/a-min(T(:))/b,13)
      plotInBox([xInl([0 1],v);yInl([0 1],v);zInl([0 1],v)],min(S(:)),max(S(:)),min(T(:)),max(T(:)));
    end
  case 'EH'
    P=subs('[a*(cos(v)+u*sin(v));b*(sin(v)-u*cos(v));c*u]',{a b},abc(1:2));
    xInl=inline(vectorize(char(P(1))),'u','v');
    yInl=inline(vectorize(char(P(2))),'u','v');
    zInl=inline(vectorize(char(P(3))),'u','v','c');
    a=abc{1};
    b=abc{2};
    c=abc{3};
    [U,V]=meshgrid([min(Z(:)) max(Z(:))]/c,linspace(0,2*pi,25));
    X=xInl(U,V);
    Y=yInl(U,V);
    Z=zInl(U,V,c);
    plot3(X',Y',Z','color',0.9*[1 1 1],'linewidth',3.0);
    Z=zInl(U,V,-c);
    plot3(X',Y',Z','color',0.9*[1 1 1],'linewidth',3.0);
  case('Eb')
  otherwise
    if strcmp(opt,'none')
      sG=linspace(min(S(:)),max(S(:)),13);
      tG=linspace(min(T(:)),max(T(:)),13);
    else
      sG=linspace(min(S(:)),max(S(:)),25);
      tG=linspace(min(T(:)),max(T(:)),25);
    end
    [S,T]=meshgrid(sG,tRg);
    X=xInl(S,T);
    Y=yInl(S,T);
    Z=zInl(S,T);
    if max(size(X))==1
      X=repmat(X,size(T));
    end
    if max(size(Y))==1
      Y=repmat(Y,size(T));
    end
    if max(size(Z))==1
      Z=repmat(Z,size(T));
    end
    if strcmp(opt,'none')
      plot3(X,Y,Z,'color',0.9*[1 1 1],'linewidth',3.0,'erase','none');
    else
      plot3(X,Y,Z,'color',0.9*[1 1 1],'linewidth',3.0);
    end
    [S,T]=meshgrid(sRg,tG);
    X=xInl(S,T);
    Y=yInl(S,T);
    Z=zInl(S,T);
    if max(size(X))==1
      X=repmat(X,size(T));
    end
    if max(size(Y))==1
      Y=repmat(Y,size(T));
    end
    if max(size(Z))==1
      Z=repmat(Z,size(T));
    end
    plot3(X',Y',Z','color',0.9*[1 1 1],'linewidth',3.0);
end

function plotInBox(P,xmin,xmax,ymin,ymax)
P0=P(:,1);
v=P(:,2)-P(:,1);
A=P0+((ymin-P0(2))/v(2))*v;
if A(1)>xmin & A(1)<xmax
  Q=A;
else
  Q=zeros(3,0);
end
A=P0+((ymax-P0(2))/v(2))*v;
if A(1)>xmin & A(1)<xmax
  Q=[Q A];
end
A=P0+((xmin-P0(1))/v(1))*v;
if A(2)>ymin & A(2)<ymax
  Q=[Q A];
end
A=P0+((xmax-P0(1))/v(1))*v;
if A(2)>ymin & A(2)<ymax
  Q=[Q A];
end
plot3(Q(1,:),Q(2,:),Q(3,:),'color',0.9*[1 1 1],'linewidth',4.0);
