function RechtStumpf3eckHM

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

clc;
close('all')
global ax
ax=[inf -inf inf -inf];
deg=pi/180;
phi=[0:360]*deg;

A=[0;0];
c=[50;0];
B=A+c;
alpha=atan(4/5)/deg;
C=SchnittGerGer(A,A+[cos(alpha);sin(alpha)],B,B+[cos(pi/2+alpha);sin(pi/2+alpha)]);
b=C-A;

Polygon([A B C A],'k','linewidth',1.5);
hold('on');
Polygon([Lotfuss(A,B,C) C],'k');
Polygon([Lotfuss(B,C,A) A],'k');
Polygon([Lotfuss(C,A,B) B],'k');

A(1)=A(1)+1.25*c(1);
B=A+c;
C=A+b;
Polygon([A B C A],'k','linewidth',1.5);
Polygon([M(A,B,C) (A+B)/2],'k');
Polygon([M(B,C,A) (B+C)/2],'k');
Polygon([M(C,A,B) (C+A)/2],'k');
Kreisbogen(B,SchnittGerGer(M(A,B,C),(A+B)/2,M(B,C,A),(B+C)/2),A,30*deg,'k');

A(1)=A(1)+1.25*c(1);
B=A+c;
alpha=30*deg;
beta=40*deg;
C=SchnittGerGer(A,A+[cos(alpha);sin(alpha)],B,B+[cos(pi-beta);sin(pi-beta)]);
b=C-A;

Polygon([A B C A],'k','linewidth',1.5);
a=C-B;
H=SchnittGerGer(Lotfuss(A,B,C),Lotfuss(A,B,C)+[-c(2);c(1)],Lotfuss(B,C,A),Lotfuss(B,C,A)+[-a(2);a(1)]);
Polygon([Lotfuss(A,B,C) Lotfuss(A,B,C)+1.2*(H-Lotfuss(A,B,C))],'k');
Polygon([A A+1.2*(H-A)],'k');
Polygon([B B+1.2*(H-B)],'k');
Polygon([C Lotfuss(A,C,B)],'k','linewidth',0.25)
Polygon([C Lotfuss(B,C,A)],'k','linewidth',0.25)

A(1)=A(1)+1.25*c(1);
B=A+c;
C=A+b;
Polygon([A B C A],'k','linewidth',1.5);
UM=SchnittGerGer(M(A,B,C),(A+B)/2,M(B,C,A),(B+C)/2);
Polygon([UM+0.6*(UM-(A+B)/2) (A+B)/2],'k');
Polygon([UM+0.2*(UM-(B+C)/2) (B+C)/2],'k');
Polygon([UM+0.2*(UM-(C+A)/2) (C+A)/2],'k');
Kreisbogen(B,UM,A,30*deg,'k');

hold('off')
axis('equal')
axis(ax)
axis('off')
print('-deps','RechtStumpf3eckHM');

function H=Lotfuss(A,B,C);
syms('t','real');
v=B-A;
H=A+double(solve(v'*(A+t*v-C),t))*(B-A);

function H=M(A,B,C);
syms('t','real');
v=B-A;
t1=double(solve(v'*(B+t*(C-B)-(A+B)/2),t));
if (t1>=0) & (t1<=1)
  H=B+t1*(C-B);
else
  H=A+double(solve(v'*(A+t*(C-A)-(A+B)/2),t))*(C-A);
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 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 Kreisbogen(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));
Polygon(repmat(M,size(t))+r*[cos(t);sin(t)],varargin{:});