%% Inkugel eines Tetraeders

clear,clc,clf
%
%% Punkte des Tetraeders
%
a = [ 0 0 0 ]; b = [ 1 0 0 ]; c = [ 0 1 0 ]; d = [ 0 0 1 ];
% a = [  5 -1 2 ]; b = [ 1 3 4 ]; c = [ -1 4 2 ]; d = [ 1 0 1 ];

%
%% Normalenvektoren der 4 Tetraederflächen
%   Rechtsregel beachten (Normalenvektoren sollen nach Innen zeigen)
%   Kann zum Bsp durch Einsetzen des Tetraederschwerpunktes S
%   in die Formel d = N x (S -X) erzielt werden. 
%   Falls d<0, dann Normalenvektor ist N = -N  notwendig.
%
n1 = cross( b-a, c-a );   % und
n1 = n1/norm(n1);         %     normalisieren

n2 = cross( d-a, b-a );   % und
n2 = n2/norm(n2);         %     normalisieren

n3 = cross( c-a, d-a );   % und
n3 = n3/norm(n3);         %     normalisieren

n4 = cross( b-c, d-c );   % und
n4 = n4/norm(n4);         %     normalisieren

%
%% Definition of Matrix A and rechter Seite f des Gleichungssystems 
% zur Berechnung des Inkreismittelpunktes M und dessen Radius r
%   Lösungsvektor u = [ m_x m_y m_z r]'
%

A = [ n1(1)  n1(2)  n1(3)  -1;
      n2(1)  n2(2)  n2(3)  -1;
      n3(1)  n3(2)  n3(3)  -1;
      n4(1)  n4(2)  n4(3)  -1; ];

f = [ dot(n1,a);  dot(n2,a);  dot(n3,a);  dot(n4,d) ];

%% und das Gleichungsystem lösen

u = A\f;

%
%% Ergebnisausgabe
%
disp(' ')
disp('Eckpunkte des Tetraeders')
disp(vertcat(a,b,c,d))
disp(' ')
disp('Inkreismittelpunkt:')
M = u(1:3);
disp(M)
disp(' ')
disp('Radius:')
Radius = u(4);
disp(Radius)



%% Grafik
%%  Tetraeder
%
%  Koordinaten
koord = horzcat(a', b', c', d');
tri = [ 1 2 3; 1 2 4; 1 3 4; 2 3 4];

trisurf(tri, koord(1,:), koord(2,:), koord(3,:) )
alpha(.5)   % transparency factor
hold


%% Kugel visualisieren

k = 5;
n = 2^k-1;

theta = pi*(-n:2:n)/n;
phi = (pi/2)*(-n:2:n)/n;

[ TT PP] = meshgrid(theta,phi);
XX = M(1)+Radius*cos(PP).*cos(TT);
YY = M(2)+Radius*cos(PP).*sin(TT);
ZZ = M(3)+Radius*sin(PP);

surf(XX,YY,ZZ)
colormap hsv;
axis equal
%pause