Unit PaulBourke;
{
http://www.alternatievewiskunde.nl/sunall/suna55.htm
}
INTERFACE

uses
  Graphics, ExtCtrls, Grafisch;

type
  iso = record
  x,y,z : double;
  end;
  isoquad = array[1..4] of iso;
var
  Beeld : TImage;

procedure Setup_Contours;
procedure Contour_Vierhoek(meest : integer ; q : isoquad);

IMPLEMENTATION

type
  rijtje = array[1..3] of integer;
var
  np : array[1..3,1..6] of integer;
  klas,doen : array[1..27] of integer;
  xd,yd,zd,fd : array[1..3] of double ;

procedure permutatie(L, n : integer ; var r : rijtje);
{
  Permutations of array r(n)
  --------------------------
  L  :=  ordinal number of permutation (input)
  n  :=  number of elements in: (input)
  r  :=  row containing the permutaton (output) }
var
  h : integer;
  k, p, q, s, m : integer;

begin
  h := 1 ;
  for k := 1 to n do
  begin
    r[k] := k ; h := h*k
  end;

{ Which item: }
  for k := 1 to n do
  begin
    p := n-k+1 ; h := h div p ;
    q := (((L-1) div h) mod p) + 1 ;
    if q = 1 then Continue;
  { Otherwise rotate data: }
    s := r[k+q-1] ;
    for m := (k+q-1) downto (k+1) do
      r[m] := r[m-1] ;
    r[k] := s ;
  end;
end;

function log2(n : integer) : integer;
{
  Logarithm base 2
  ---------------- }
var
  k, L : integer;
begin
  log2 := 0;
  if n <= 0 then Exit;
  for k := 1 to 32 do
  begin
    n := n div 2 ; L := k;
    if n = 0 then Break ;
  end;
  log2 := L;
end;

procedure Controle;
var
  i,j,k : integer;
begin
  for j := 1 to 6 do
  begin
    for i := 1 to 3 do
    begin
      Write(np[i,j]:3);
    end;
    Writeln;
  end;
  for k := 1 to 27 do
  begin
    Writeln(klas[k]:3,doen[k]:3);
  end;
end;

procedure Setup_Contours;
{
  Initialize np[3,6],klas[27],doen[27]
  ------------------------------------
  All possible cases to be classified:
  === ==< ==> =<= =<< =<> =>= =>< =>>
  <== <=< <=> <<= <<< <<> <>= <>< <>>
  >== >=< >=> ><= ><< ><> >>= >>< >>>
  Being equivalent to counting base 3,
  where 0 : = , 1 : < , 2 : > .
  Resulting in 6 equivalence classes:
       1   2   3   4   5   6    keer
      === ==< =<< =<> <<< <<>    1
          ==> =>> =>< >>> >><    2
          =<= <=< <=>     <><    1
          =>= >=> >=<     ><>    2
          <== <<= <>=     ><<    1
          >== >>= ><=     <>>    2
}
var
  L, i, j, n, m, keer : integer;
  num, rij : rijtje;
begin
  for L := 1 to 6 do
  begin
    permutatie(L, 3, rij) ;
    for j := 1 to 3 do
    begin
      np[j,L] := rij[j] ;
    end;
  end;

  for L := 1 to 27 do
  begin
    klas[L] := 0 ; doen[L] := 0 ;
  end;

{ Equivalence classes: }
  n := 0 ;
  for L := 1 to 27 do
  begin
    if klas[L] > 0 then Continue;
    n := n+1 ; m := L-1 ;
  { Decimals base 3: }
    for i := 1 to 3 do
    begin
      num[3-i+1] := m mod 3 ;
      m := m div 3 ;
    end;
    for keer := 1 to 2 do
    begin
      for i := 1 to 6 do
      begin
      { Inverse permutation: }
        for j := 1 to 3 do
          rij[np[j,i]] := j ;
        m := 0 ;
        for j := 1 to 3 do
          m := m*3+num[rij[j]] ;
        m := m+1 ;
      { Assign equivalent elements: }
        if klas[m] = 0 then klas[m] := n ; { Class }
        if doen[m] = 0 then doen[m] := i ; { Permutation }
      end;

      for i := 1 to 3 do
        num[i] := (3-num[i]) mod 3 ;
      m := 0 ;
      for j := 1 to 3 do
        m := m*3 + num[rij[j]] ;
      m := m+1 ; if klas[m] > 0 then Break ;
    end; { < and > exchanged 2nd time }
  end;
{ Controle; }
end;

procedure onder_boven(q : double ; most : integer ; var n1, n2 : integer);
{
  Determine lower & upper boundary values
  --------------------------------------- }
var
  m, grens, k : integer;
begin
  n1 := 0 ; n2 := most ; m := most ;
  grens := log2(most) ;
  for k := 1 to grens do
  begin
    m := m div 2 ;
    n1 := n1 + m ; n2 := n2 - m ;
    if q <= n1/most then n1 := n1 - m ;
    if q >= n2/most then n2 := n2 + m ;
  end;
end;

procedure Contour_Driehoek(zero : double);
{
  Contourlines in a Triangle
  -------------------------- }
var
  k, more, merk, keus : integer;
  np1, np2, np3 : integer;
  x1, x2, x3 : double;
  y1, y2, y3 : double;
  f1, f2, f3 : double;
  noemer, xc, yc : double ;
begin
  for k := 1 to 3 do
    fd[k] := zd[k] - zero ;
{ Coding each case: }
  more := 0 ;
  for k := 1 to 3 do
  begin
    if fd[k] = 0 then more := more*3 ;
    if fd[k] < 0 then more := more*3+1 ;
    if fd[k] > 0 then more := more*3+2 ;
  end;
  merk := klas[more+1] ;
  if (merk = 1) or (merk = 3) or (merk = 5) then Exit ;
{          ===           =<<           <<<            }

{ Permutation of (1,2,3): }
  keus := doen[more+1] ;
  np1 := np[1,keus] ; np2 := np[2,keus] ; np3 := np[3,keus] ;
  x1 := xd[np1] ; x2 := xd[np2] ; x3 := xd[np3] ;
  y1 := yd[np1] ; y2 := yd[np2] ; y3 := yd[np3] ;
  f1 := fd[np1] ; f2 := fd[np2] ; f3 := fd[np3] ;

  Case merk of

  1, 3, 5 : Exit ; { === , =<< , <<< }

  2 : { line through side: ==< }
      with Beeld.Canvas do
      begin
        MoveTo(x2i(x1),y2j(y1)) ;
        LineTo(x2i(x2),y2j(y2)) ;
        Exit;
      end;

  4 : { line through corner: =<> }
      with Beeld.Canvas do
      begin
        MoveTo(x2i(x1),y2j(y1)) ;
        noemer := 1/(f2-f3) ;
        xc := (f2*x3-f3*x2)*noemer;
        yc := (f2*y3-f3*y2)*noemer;
        LineTo(x2i(xc),y2j(yc)) ;
        Exit;
      end;

  6 : { line through two sides: <<> }
      with Beeld.Canvas do
      begin
        noemer := 1/(f1-f3) ;
        xc := (f1*x3-f3*x1)*noemer;
        yc := (f1*y3-f3*y1)*noemer;
        MoveTo(x2i(xc),y2j(yc)) ;
        noemer := 1/(f2-f3) ;
        xc := (f2*x3-f3*x2)*noemer;
        yc := (f2*y3-f3*y2)*noemer;
        LineTo(x2i(xc),y2j(yc)) ;
        Exit;
      end;
  end;
end;

procedure Contour_Vierhoek(meest : integer ; q : isoquad);
{
  Contourlines in a Quadrilateral
  ------------------------------- }
const
  volg : array[1..4] of integer = (1,2,4,3);
var
  i,j,k,n,nq1,nq2,nop : integer;
  qmin,qmax,zero : double;
begin
  qmin := 1 ; qmax := 0 ;
  for i := 1 to 4 do
  begin
    if q[i].z < qmin then qmin := q[i].z ;
    if q[i].z > qmax then qmax := q[i].z ;
  end;
  onder_boven(qmin, meest, nq1, nop) ;
  onder_boven(qmax, meest, nop, nq2) ;

{ Loop through Triangles: }
  xd[3] := 0.25*(q[1].x+q[2].x+q[3].x+q[4].x);
  yd[3] := 0.25*(q[1].y+q[2].y+q[3].y+q[4].y);
  zd[3] := 0.25*(q[1].z+q[2].z+q[3].z+q[4].z);
  for i := 1 to 4 do
  begin
    for j := 1 to 2 do
    begin
      k := volg[((i+j-2) mod 4) + 1];
      xd[j] := q[k].x ; yd[j] := q[k].y ;
      zd[j] := q[k].z ;
    end;
  { Loop through Levels: }
    for n := nq1 to nq2 do
    begin
      zero := n/meest ;
    { Elementary contours: }
      Contour_Driehoek(zero);
    end;
  end;
end;

END.