function SHLK

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

clc;
global ax
ax=[inf -inf inf -inf];
P=Konstr(ax,inline('','dummy1','dummy2'));
a=[min(P(1,:)) max(P(1,:)) min(P(2,:)) max(P(2,:))];
Konstr(a+[(a(2)-a(1))*[-0.1 0.1] (a(4)-a(3))*[-0.1 0.1]],@fig2eps);

function fig2eps(ax,Name)
axis('off')
axis('equal')
axis(ax)
print('-depsc','-loose',Name)

function P=Konstr(ax,pr)
close('all')
deg=pi/180;
H=[9;9];
L=[4;10];
K=[11;15];
alphaHL=22*deg;
alphaHK=40*deg;

Name='SHLK-Vgb';
figure('numbertitle','off','name',Name);
drawPolygon(@plot,[H L K H],'k','linewidth',1.5);
Schrift(H,'\itH','FontSize',40,'HorizontalAlignment','center','VerticalAlignment','top');
Schrift(L,'\itL','FontSize',40,'HorizontalAlignment','center','VerticalAlignment','top');
Schrift(K,'\itK','FontSize',40,'HorizontalAlignment','center','VerticalAlignment','bottom');
feval(pr,ax,Name);

Name='SHLK-HK';
figure('numbertitle','off','name',Name);
drawPolygon(@plot,[H L K H],'k','linewidth',1.5);
Schrift(H,'\itH','FontSize',27,'HorizontalAlignment','center','VerticalAlignment','top');
Schrift(L,'\itL','FontSize',27,'HorizontalAlignment','center','VerticalAlignment','top');
Schrift(K,'\itK','FontSize',27,'HorizontalAlignment','center','VerticalAlignment','bottom');
hold('on');
KrHL=OK(H,K,alphaHK,[0 1 0]);
feval(pr,ax,Name);

Name='SHLK';
figure('numbertitle','off','name',Name);
drawPolygon(@plot,[H L K H],'k','linewidth',1.5);
Schrift(H,'\itH','FontSize',13,'HorizontalAlignment','center','VerticalAlignment','top');
Schrift(L,'\itL','FontSize',13,'HorizontalAlignment','center','VerticalAlignment','top');
Schrift(K,'\itK','FontSize',13,'HorizontalAlignment','center','VerticalAlignment','bottom');
hold('on');
KrHL=OK(H,L,alphaHL,[1 0 1]);
feval(pr,ax,[Name '-HL']);

KrHK=OK(H,K,alphaHK,[0 1 0]);
S=SchnittOK(KrHL,KrHK);

alphaLK=alphaHL+alphaHK
for s=1:size(S,2);
  drawPolygon(@plot,[L S(:,s) K],'--k');
  a=L-S(:,s);
  b=K-S(:,s);
  if abs(dot(a/norm(a),b/norm(b))-cos(alphaLK))<sqrt(eps)
    drawPolygon(@plot,S(:,s),'.','color',0.00*[1 1 1],'MarkerSize',30)
  else
    drawPolygon(@plot,S(:,s),'.','color',0.75*[1 1 1],'MarkerSize',30)
  end
end

P=[S H L K KrHL(1).M KrHL(2).M KrHK(1).M KrHK(2).M];
hold('off')
feval(pr,ax,Name);

function S=SchnittOK(K1,K2)
co=[cos(K1(1).alpha) cos(K2(1).alpha)];
A=[K1(1).A K2(1).A];
B=[K1(1).B K2(1).B];
S=[];
delta=sqrt(eps);
for a=1:2
  for b=1:2
    [T(:,1),T(:,2)]=SchnittKreisKreis(K1(a),K2(b));
    for c=1:2
      vA=A-T(:,[c c]);
      vB=B-T(:,[c c]);
      nA=sqrt(dot(vA,vA));
      nB=sqrt(dot(vB,vB));
      if all([nA nB])>delta
        vA=vA/diag(nA);
        vB=vB/diag(nB);
        if all(abs(dot(vA,vB)-co)<delta)
          S=[S T(:,c)];
        end
      end
    end
  end
end

function K=OK(A,B,alpha,col)
deg=pi/180;
drawPolygon(@plot,[A B],'k','linewidth',1.5);
t0=cart2pol(B(1)-A(1),B(2)-A(2));
tA=t0+90*deg-alpha;
tB=t0+90*deg+alpha;
K(1).M=SchnittGerGer(A,A+[cos(tA);sin(tA)],B,B+[cos(tB);sin(tB)]);
K(1).r=norm(K(1).M-A);
K(1).A=A;
K(1).B=B;
K(1).alpha=alpha;
K(2)=K(1);
K(2).M=SchnittGerGer(A,A+[cos(tB);sin(tB)],B,B+[cos(tA);sin(tA)]);
drawPolygon(@plot,[A A+1.1*(K(1).M-A)],'color',col)
drawPolygon(@plot,[B B+1.1*(K(1).M-B)],'color',col)
drawPolygon(@plot,[A A+1.1*(K(2).M-A)],'color',col)
drawPolygon(@plot,[B B+1.1*(K(2).M-B)],'color',col)
v=(K(1).M-K(2).M)/2;
drawCircle(@plot,B,K(1).M,A,0,'color',col,'linewidth',1.5);
drawCircle(@plot,A,K(1).M,B,0,'color',col);
drawCircle(@plot,A,K(2).M,B,0,'color',col,'linewidth',1.5);
drawCircle(@plot,B,K(2).M,A,0,'color',col);

function S=SchnittGerGer(A,B,C,D)
syms('t','real');
n=[C(2)-D(2);D(1)-C(1)];
S=A+double(solve(n'*(A+t*(B-A)-C),t))*(B-A);

function [S1,S2]=SchnittKreisKreis(K1,K2)
syms('t','real');
P=K2.M+K2.r*[cos(t);sin(t)];
tSol=double(solve((K1.M-P)'*(K1.M-P)-K1.r^2,t));
S1=K2.M+K2.r*[cos(tSol(1));sin(tSol(1))];
if abs(imag(S1))<sqrt(eps)*(1+abs(S1))
  S1=real(S1);
end
S2=K2.M+K2.r*[cos(tSol(2));sin(tSol(2))];
if abs(imag(S2))<sqrt(eps)*(1+abs(S2))
  S2=real(S2);
end

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 C=drawCircle(draw,A,M,B,offset,varargin)
deg=pi/180;
[t1,r1]=cart2pol(A(1)-M(1),A(2)-M(2));
[t2,r2]=cart2pol(B(1)-M(1),B(2)-M(2));
if abs(r1-r2)>sqrt(eps)
  error('Kreisbogen: A und B liegen nicht auf einem Kreis mit Mittelpunkt M')
end
r=(r1+r2)/2;
t1=t1-offset;
t2=t2+offset;
if t2<t1
  t2=t2+2*pi;
end
t=linspace(t1,t2,round((t2-t1)/deg));
C=repmat(M,size(t))+r*[cos(t);sin(t)];
drawPolygon(draw,C,varargin{:});

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

function r=Rot(a,b)
deg=pi/180;
r=atan2(b(2)-a(2),b(1)-a(1))/deg;
if r<0
  r=r+180;
end
if r>90
  r=r-180;
end
