function ESV

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

close('all');
clc
deg=pi/180;
OrthAxo(25*deg,35*deg,40*deg,'ESV-Tri',1);
OrthAxo(asin(1/8),acos(1/4*sqrt(7))-asin(1/3*sqrt(2)),asin(1/4*sqrt(7))+asin(1/3*sqrt(2)),'ESV-Di1',[]);
OrthAxo(asin(1/4*sqrt(7)),asin(1/4*sqrt(7))+asin(1/3*sqrt(2)),45*deg,'ESV-Di2',1);
OrthAxo(30*deg,45*deg,45*deg,'ESV-Iso',1);

function OrthAxo(alpha,gamma1,gamma2,Name,spez)
global ax
figure('numbertitle','off','name',Name);
ax=[inf -inf inf -inf];
axo.alpha=alpha;
axo.p1=20;
axo.p2=29;
axo.p3=15;
axo.xy=Riss(gamma1,15,17,20);
axo.yz=Riss(gamma2,10,20,15);
hold('on')
axo=initAchsen(axo)
plotGA(axo,[12 12 12 12   12 12
             0 15 15  7.5  0  0
             0  0  7 14    7  0]);
plotGA(axo,[ 0  0  0  
            15 15  7.5
             0  7 14  ]);
for z=[0 7]
  plotGA(axo,[ 0 12
              15 15
               z z]);
end        
plotGA(axo,[ 0   12
             7.5  7.5
            14   14]);
if ismember(1,spez)
  plotGA(axo,[ 0   0 12
               7.5 0  0
              14   7  7]);
else
  [G,A,S]=GAS(axo,[ 0   0 12
                    7.5 0  0
                   14   7  7]);
  drawPolygon(@plot,G,'k','linewidth',2.0)
  drawPolygon(@plot,A,'k','linewidth',2.0)
end            
hold('off')
axis('off')
axis('equal')
axis(ax)
print('-deps',Name)

function R=Riss(gamma,dO,p1,p2)
R.gamma=gamma;
R.dO=dO; 
R.p1=p1; 
R.p2=p2;

function Axo=initAchsen(axo)
deg=pi/180;
Axo=axo;
Axo.x=[-cos(axo.alpha);-sin(axo.alpha)];
Axo.z=[0;1];
L=max([axo.p1;axo.p2;axo.p3])/10;
[Axo.xy.O,Axo.xy.x,Axo.xy.y]=plotAR(pi/2,pi+axo.xy.gamma,axo.xy,L,'x''','y''');
[Axo.yz.O,Axo.yz.y,Axo.yz.z]=plotAR(pi+axo.alpha,1.5*pi+axo.alpha+axo.yz.gamma,axo.yz,L,'y''''','z''''');
[G,A,S]=GAS(Axo,[0;1;0]);
Axo.y=S/norm(S);
O=[0;0];
P1=Axo.x*axo.p1;
P2=Axo.y*axo.p2;
P3=Axo.z*axo.p3;
drawPolygon(@plot,[O P1],'k');
drawPolygon(@plot,[O P2],'k');
drawPolygon(@plot,[O P3],'k');
Spitze(P1,P1,L,30*deg);
Spitze(P2,P2,L,30*deg);
Spitze(P3,P3,L,30*deg);
Schrift(P1,'\itx ','HorizontalAlignment','right','VerticalAlignment','middle')
Schrift(P2,'\ity','HorizontalAlignment','left','VerticalAlignment','middle')
Schrift(P3,'\itz','HorizontalAlignment','center','VerticalAlignment','bottom')

function [O1,e1,e2]=plotAR(w,w1,R,L,s1,s2);
deg=pi/180;
s=sign(w1);
w1=mod(w1,2*pi);
O=[0;0];
O1=-[cos(w);sin(w)]*R.dO;
w2=mod(w1+0.5*s*pi,2*pi);
e1=[cos(w1);sin(w1)];
e2=[cos(w2);sin(w2)];
P1=O1+e1*R.p1;
P2=O1+e2*R.p2;
drawPolygon(@plot,[O O1],':k')
drawPolygon(@plot,[P2 O1 P1],'k')
Spitze(P1,P1-O1,L,30*deg);
Spitze(P2,P2-O1,L,30*deg);
if w1>1.5*pi | w1<0.5*pi
  Schrift(P1,['\it' s1],'HorizontalAlignment','left','VerticalAlignment','middle')
else
  Schrift(P1,['\it' s1],'HorizontalAlignment','right','VerticalAlignment','middle')
end
if w2<pi
  Schrift(P2,['\it' s2],'HorizontalAlignment','center','VerticalAlignment','bottom')
else
  Schrift(P2,['\it' s2],'HorizontalAlignment','center','VerticalAlignment','top')
end

function [G,A,S]=GAS(axo,P) 
G=ProjRiss(axo.xy.O,axo.xy.x,axo.xy.y,P(1:2,:));
A=ProjRiss(axo.yz.O,axo.yz.y,axo.yz.z,P(2:3,:));
n=[-axo.x(2);axo.x(1)];
S=zeros(2,size(P,2));
for k=1:size(P,2)
  S(:,k)=G(:,k)-(((G(:,k)-A(:,k))'*n)/(axo.z'*n))*axo.z;
end

function [G,A,S]=plotGA(axo,P)
[G,A,S]=GAS(axo,P);
for k=1:size(P,2)
  drawPolygon(@plot,[G(:,k) S(:,k) A(:,k)],'--k','linewidth',0.25)
end
drawPolygon(@plot,G,'k','linewidth',2.0)
drawPolygon(@plot,A,'k','linewidth',2.0)
drawPolygon(@plot,S,'k','linewidth',2.0)

function Q=ProjRiss(o,e1,e2,P)
v1=repmat(P(1,:),2,1).*repmat(e1,1,size(P,2));
v2=repmat(P(2,:),2,1).*repmat(e2,1,size(P,2));
O=repmat(o,1,size(P,2));
Q=O+v1+v2;

function drawPolygon(draw,P,varargin);
global ax
feval(draw,P(1,:),P(2,:),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

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

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