Unit Algemeen;
{
  GENERAL ROUTINES
}
INTERFACE

Uses SysUtils;

type
  punt = record
    x,y : double;
  end;
  punten = array of punt;
  drie = record
    A,B,C : integer;
  end;
  drietal = array of drie;
  integers = array of integer;

procedure TijdMeting(afdruk : string); { Computing Time Measurement }
function determinant(a,b,c : punt) : double; { Twice area of triangle }
function binnen(a,b,c,p : punt) : boolean; { (p) is inside triangle (a,b,c) }
procedure opnieuw(N : integer; var rij : punten); { New Set of Random Sites }
procedure omcirkel(A,B,C : punt; var m : punt); { Circumscribed Circle with midpoint (m) }

IMPLEMENTATION

var
  eerder : integer;

procedure TijdMeting(afdruk : string);
{
  StopWatch
}
var
  heden : integer;
begin
  if afdruk = '' then
  begin
    eerder := DateTimeToTimeStamp(Now).Time;
    Exit;
  end;
  heden := DateTimeToTimeStamp(Now).Time;
  Writeln(afdruk+' = '+IntToStr(heden-eerder)+' ms');
end;

function determinant(a,b,c : punt) : double;
{
  (Two times) Area of Triangle
}
begin
  determinant := (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
end;

function min(a,b,c : double) : double;
{
  Minimum of (a,b,c)
}
var
  m : double;
begin
  m := a;
  if b < a then m := b;
  if c < m then m := c;
  min := m;
end;

function binnen(a,b,c,p : punt) : boolean;
{
  (p) is inside triangle (a,b,c)
}
var
  det,xi,eta : double;
begin
  det := (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
  xi := ((c.y-a.y)*(p.x-a.x)-(c.x-a.x)*(p.y-a.y))/det;
  eta := ((b.x-a.x)*(p.y-a.y)-(b.y-a.y)*(p.x-a.x))/det;
  binnen := min(xi,eta,1-xi-eta) > 0;
end;

procedure opnieuw(N : integer; var rij : punten);
{
  Random Sites (rij)
}
var
  k : integer;
begin
{ Random points in the plane }
  SetLength(rij,N);
  for k := 0 to N-1 do
  begin
    rij[k].x := Random; rij[k].y := Random;
  end;
end;

procedure omcirkel(A,B,C : punt; var m : punt);
{
  Cirumscribed circle of triangle (A,B,C)
  output its (m) = midpoint, (R) = radius
}
var
  x0,x1,x2,y0,y1,y2 : double;
  rhs1,rhs2,det,p,q : double;
begin
  x0 := A.x; y0 := A.y;
  x1 := B.x; y1 := B.y;
  rhs1 := sqr(x1)-sqr(x0)+sqr(y1)-sqr(y0);
  x2 := C.x; y2 := C.y;
  rhs2 := sqr(x2)-sqr(x0)+sqr(y2)-sqr(y0);
  det := (x1-x0)*(y2-y0)-(y1-y0)*(x2-x0);
  p := ((y2-y0)*rhs1-(y1-y0)*rhs2)/(2*det);
  q := (-(x2-x0)*rhs1+(x1-x0)*rhs2)/(2*det);
  m.x := p; m.y := q;
end;

END.