function StrahlensatzB

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

clc;
close('all')
global ax
ax=[inf -inf inf -inf];
deg=pi/180;
t1=220*deg;
t2=t1+120*deg;
v1=[cos(t1);sin(t1)]
v2=[cos(t2);sin(t2)]
a1=5; a2=7.5; b1=2;b2=a2/a1*b1;
S=[0;0];
A1=S+a1*v1; B1=S-b1*v1; A2=S+a2*v2; B2=S-b2*v2;
Polygon([A1-0.1*(B1-A1) A1+1.1*(B1-A1)],'k')
hold('on')

Polygon([A2-0.1*(B2-A2) A2+1.1*(B2-A2)],'k')
v=A2-A1;
Polygon([A1-0.1*v A2+0.1*v],'k')
Polygon([B1+0.1*v B2-0.1*v],'k')

J=inline('[-v(2);v(1)]','v');
q=0.2;
n=J(v)/norm(v);
A1inn=SchnittGerGer(A1+q*n,A2+q*n,S+q*J(v1),A1+q*J(v1));
A2inn=SchnittGerGer(A1+q*n,A2+q*n,S-q*J(v2),A2-q*J(v2));
Sinn=SchnittGerGer(S+q*J(v1),A1+q*J(v1),S-q*J(v2),A2-q*J(v2));
PolygonM([Sinn A1inn A2inn Sinn],q,'.','color',0.5*[1 1 1],'MarkerSize',15)
B1aus=SchnittGerGer(B1-q*n,B2-q*n,S-q*J(v1),A1-q*J(v1));
B2aus=SchnittGerGer(B1-q*n,B2-q*n,S+q*J(v2),A2+q*J(v2));
Saus=SchnittGerGer(S-q*J(v1),A1-q*J(v1),S+q*J(v2),A2+q*J(v2));
PolygonM([Saus B1aus B2aus Saus],q,'o','color',0.5*[1 1 1])

Schrift((S+A1)/2,'\ita_1','HorizontalAlignment','center','VerticalAlignment','bottom','rotation',Rot(A1,S))
Schrift((S+B1)/2,'\itc_1','HorizontalAlignment','center','VerticalAlignment','top','rotation',Rot(A1,S))
Schrift((S+A2)/2,'\ita_2','HorizontalAlignment','center','VerticalAlignment','bottom','rotation',Rot(A2,S))
Schrift((S+B2)/2,'\itc_2','HorizontalAlignment','center','VerticalAlignment','top','rotation',Rot(A2,S))
Schrift((A1+A2)/2,'\ite','HorizontalAlignment','center','VerticalAlignment','top','rotation',Rot(A1,A2))
Schrift((B1+B2)/2,'\itf','HorizontalAlignment','center','VerticalAlignment','bottom','rotation',Rot(A1,A2))

hold('off')
axis('equal')
axis(ax)
axis('off')
print('-depsc','StrahlensatzB');

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 Polygon(P,varargin);
global ax
plot(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 PolygonM(P,d,varargin)
for k=2:size(P,2)
  v=P(:,k)-P(:,k-1);
  t=linspace(0,1,round(norm(v)/d));
  Polygon(repmat(P(:,k-1),1,length(t))+repmat(t,length(P(:,k-1)),1).*repmat(v,1,length(t)),varargin{:})
end

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 Schrift(P,txt,varargin)
text(P(1),P(2),txt,'FontName','TimesNewRoman','FontSize',20,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
