Unit IsoSurf;
{
        program Testen;

        Uses IsoSurf;

        begin
          Bol;
          Kubus;
        end.
}
INTERFACE

type
  row8d = array[1..8] of double;
  row6d = array[1..6] of double;

var
  uit : text;

procedure Setup_Contours;
procedure Contour_Blok(zero : double ; xq,yq,zq,fq : row8d);
procedure Contour_Prisma(zero : double ; xp,yp,zp,fp : row6d);
procedure Bol ; procedure Kubus ;

IMPLEMENTATION

type
  rijtje = array[1..4] of byte;
  row4d = array[1..8] of double;

const { Even permutations only: }
  np : array[1..12,1..4] of byte =
       ((1,2,3,4),(2,1,4,3),(3,4,1,2),(4,3,2,1),
        (2,3,1,4),(1,4,2,3),(4,1,3,2),(3,2,4,1),
        (3,1,2,4),(4,2,1,3),(1,3,4,2),(2,4,3,1));
var
  klas, move : array[1..81] of byte;

{************************************}

function Volume(p,q,r : row4d) : double;
var
  a,b,c : array[1..3] of double;
  k : byte;
begin
  for k := 1 to 3 do
  begin
    a[k] := p[k+1]-p[1];
    b[k] := q[k+1]-q[1];
    c[k] := r[k+1]-r[1];
  end;
  volume := a[1]*b[2]*c[3] - a[1]*b[3]*c[2]
          - a[2]*b[1]*c[3] + a[2]*b[3]*c[1]
          + a[3]*b[1]*c[2] - a[3]*b[2]*c[1];
end;

{************************************}

procedure Setup_Contours;
{
  Initialize klas[81],move[81]
  ----------------------------
  All possible cases to be classified
  being equivalent to counting base 3,
  where 0 : = , 1 : < , 2 : > .
  Resulting in fifteen Equivalence Classes:
    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
  ==== ===< ===> ==<< ==<> ==>> =<<< =<<> =<>> =>>> <<<< <<<> <<>> <>>> >>>>
}
var
  L, i, j, n, m, keer : byte;
  num, rij : rijtje;

begin
  for L := 1 to 81 do
  begin
    klas[L] := 0 ; move[L] := 0 ;
  end;

{ Equivalence classes: }
  n := 0 ;
  for L := 1 to 81 do
  begin
  { Already covered: }
    if klas[L] > 0 then Continue;
    n := n + 1 ; m := L - 1 ;
  { Decimals base 3: }
    for i := 1 to 4 do
    begin
      num[4-i+1] := m mod 3 ;
      m := m div 3 ;
    end;

    for i := 1 to 12 do
    begin
    { Inverse permutation: }
      for j := 1 to 4 do
        rij[np[i,j]] := j ;

      m := 0 ;
      for j := 1 to 4 do
        m := m*3 + num[rij[j]] ;
      m := m + 1 ;

    { Assign equivalent elements: }
      if klas[m] = 0 then klas[m] := n ; { Class }
      if move[m] = 0 then move[m] := i ; { Permutation }
    end;

  end;
end;

{************************************}

procedure Contour_Tetraeder(x,y,z,f : row4d);
{
  Contours in a Tetrahedron
  ------------------------- }
var
{ Permuted & calculated: }
  cv : array[1..3,1..4] of double ;
  ca,cb : array[1..3,1..3] of double ;
  fv : array[1..4] of double ;
  i,k, meer, merk, genoeg : byte;
  fac : double;

function afstand(L : byte) : double;
var k : byte ;
var h : double;
begin
  h := 0;
  for k := 1 to 3 do
    h := h + sqr(ca[k,L]-cb[k,L]);
  afstand := sqrt(h);
end;

begin
{ Coding each case: }
  meer := 0 ;
  for k := 1 to 4 do
  begin
    if f[k] = 0 then meer := meer*3 ;
    if f[k] < 0 then meer := meer*3+1 ;
    if f[k] > 0 then meer := meer*3+2 ;
  end;

  merk := klas[meer+1];
  if (merk < 1) or (merk > 15) then
  begin
  { Should Never happen: }
    Writeln(uit,'NO contour_tetraeder');
    Close(uit);
    Halt;
  end;

{ Permutation of case: }
  for i := 1 to 4 do
  begin
    k := np[move[meer+1],i] ;
    cv[1,i] := x[k] ;
    cv[2,i] := y[k] ;
    cv[3,i] := z[k] ;
    fv[i] := f[k] ;
  end;

  genoeg := 1;
{
    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
  ==== ===< ===> ==<< ==<> ==>> =<<< =<<> =<>> =>>> <<<< <<<> <<>> <>>> >>>>
}
  Case merk of

  1, 4, 6, 7, 10, 11, 15 : Exit ; { ==== ==<< ==>> =<<< =>>> <<<< >>>> }

  2 : begin { coplanar with boundary: ===< }

      for i := 1 to 3 do
      for k := 1 to 3 do
        ca[k,4-i] := cv[k,i];

      end;

  3 : begin { coplanar with boundary: ===> }

      for i := 1 to 3 do
      for k := 1 to 3 do
        ca[k,i] := cv[k,i];

      end;

  5 : begin { through one edge: ==<> }

      for i := 1 to 2 do
      for k := 1 to 3 do
        ca[k,i] := cv[k,i];
      fac := - fv[3]/(fv[4]-fv[3]) ;
      for k := 1 to 3 do
        ca[k,3] := cv[k,3] + fac*(cv[k,4]-cv[k,3]);

      end;

  8 : begin { through one vertex: =<<> }

      for k := 1 to 3 do
        ca[k,1] := cv[k,1];
      for i := 2 to 3 do
      begin
        fac := - fv[i]/(fv[4]-fv[i]);
        for k := 1 to 3 do
          ca[k,i] := cv[k,i] + fac*(cv[k,4]-cv[k,i]);
      end;

      end;

  9 : begin { through one vertex: =<>> }

      for k := 1 to 3 do
        ca[k,1] := cv[k,1];
      for i := 3 to 4 do
      begin
        fac := - fv[i]/(fv[2]-fv[i]);
        for k := 1 to 3 do
          ca[k,6-i] := cv[k,i] + fac*(cv[k,2]-cv[k,i]);
      end;

      end;

 12 : begin { triangular surface: <<<> }

      for i := 1 to 3 do
      begin
        fac := - fv[i]/(fv[4]-fv[i]) ;
        for k := 1 to 3 do
          ca[k,i] := cv[k,i] + fac*(cv[k,4]-cv[k,i]) ;
      end;

      end;

 14 : begin { triangular surface: <>>> }

      for i := 2 to 4 do
      begin
        fac := - fv[i]/(fv[1]-fv[i]) ;
        for k := 1 to 3 do
          ca[k,i-1] := cv[k,i] + fac*(cv[k,1]-cv[k,i]) ;
      end;

      end;

 13 : begin { quadrilateral surface: <<>> }

      for i := 1 to 2 do
      begin
        fac := - fv[i]/(fv[4]-fv[i]) ;
        for k := 1 to 3 do
          ca[k,i] := cv[k,i] + fac*(cv[k,4]-cv[k,i]) ;
      end ;
      for i := 1 to 2 do
      begin
        fac := - fv[i]/(fv[3]-fv[i]) ;
        for k := 1 to 3 do
          cb[k,3-i] := cv[k,i] + fac*(cv[k,3]-cv[k,i]) ;
      end;

      { giving two triangles } genoeg := 2;

      if afstand(1) < afstand(2) then
      begin
        for k := 1 to 3 do
        begin
          ca[k,3] := cb[k,1];
          cb[k,3] := ca[k,1];
        end;
      end else begin
        for k := 1 to 3 do
        begin
          ca[k,3] := cb[k,2];
          cb[k,3] := ca[k,2];
        end;
      end;

      end;

  end;

  for i := 1 to 3 do
  begin
    for k := 1 to 3 do
      Write(uit,ca[k,4-i],' ');
    Writeln(uit);
  end;
  Writeln(uit);

  if genoeg = 1 then Exit;

  for i := 1 to 3 do
  begin
    for k := 1 to 3 do
      Write(uit,cb[k,4-i],' ');
    Writeln(uit);
  end;
  Writeln(uit);

end;

{************************************}

procedure Contour_Blok;

var
  i : byte;
  xt,yt,zt,ft : row4d;
  qmin, qmax : double;

procedure Voor_Tetraeder(welke : byte ; vq : row8d ; var vt : row4d);

const
  r : array[1..6,1..4] of byte =
  ((1,3,7,5),(2,6,8,4),(1,5,6,2),(3,4,8,7),(1,2,4,3),(7,8,6,5));

var
  a,b,c : byte;

begin
  vt[4] := ( vq[1] + vq[2] + vq[3] + vq[4]
           + vq[5] + vq[6] + vq[7] + vq[8] ) / 8 ;

  a := (welke-1) div 4 + 1;
  vt[1] := ( vq[r[a,1]] + vq[r[a,2]] + vq[r[a,3]] + vq[r[a,4]] ) / 4 ;

  b := (welke-1) mod 4 + 1;
  c := b mod 4 + 1;
  vt[2] := vq[r[a,b]] ; vt[3] := vq[r[a,c]] ;
end;

begin
  qmin := fq[1] ; qmax := fq[1];
  for i := 2 to 8 do
  begin
    if fq[i] < qmin then qmin := fq[i] ;
    if fq[i] > qmax then qmax := fq[i] ;
  end;
  if (zero < qmin) or (zero > qmax) then Exit;
  for i := 1 to 8 do
    fq[i] := fq[i] - zero;

{ Loop through Tetrahedra: }

  for i := 1 to 24 do
  begin
    Voor_Tetraeder(i,xq,xt);
    Voor_Tetraeder(i,yq,yt);
    Voor_Tetraeder(i,zq,zt);
    Voor_Tetraeder(i,fq,ft);

  { Elementary contours: }
    Contour_Tetraeder(xt,yt,zt,ft);
  end;

end;

{************************************}

procedure Contour_Prisma;

var
  i : byte;
  xt,yt,zt,ft : row4d;
  pmin, pmax : double;

procedure Voor_Tetraeder(welke : byte ; vp : row6d ; var vt : row4d);

const
  r : array[1..3,1..4] of byte =
  ((2,1,4,5),(3,2,5,6),(1,3,6,4));

var
  a,b,c : byte;

begin
  vt[4] := ( vp[1] + vp[2] + vp[3]
           + vp[4] + vp[5] + vp[6] ) / 6 ;

  if welke <= 12 then
  begin
    a := (welke-1) div 4 + 1;
    vt[1] := ( vp[r[a,1]] + vp[r[a,2]] + vp[r[a,3]] + vp[r[a,4]] ) / 4 ;

    b := (welke-1) mod 4 + 1;
    c := b mod 4 + 1;
    vt[2] := vp[r[a,b]] ; vt[3] := vp[r[a,c]] ;
    Exit;
  end;

  if welke = 13 then
  begin
    vt[1] := 1 ; vt[2] := 2 ; vt[3] := 3 ;
  end;
  if welke = 14 then
  begin
    vt[1] := 6 ; vt[2] := 5 ; vt[3] := 4 ;
  end;

end;

begin
  pmin := fp[1] ; pmax := fp[1];
  for i := 2 to 6 do
  begin
    if fp[i] < pmin then pmin := fp[i] ;
    if fp[i] > pmax then pmax := fp[i] ;
  end;
  if (zero < pmin) or (zero > pmax) then Exit;
  for i := 1 to 6 do
    fp[i] := fp[i] - zero;

{ Loop through Tetrahedra: }

  for i := 1 to 14 do
  begin
    Voor_Tetraeder(i,xp,xt);
    Voor_Tetraeder(i,yp,yt);
    Voor_Tetraeder(i,zp,zt);
    Voor_Tetraeder(i,fp,ft);

  { Elementary contours: }
    Contour_Tetraeder(xt,yt,zt,ft);
  end;

end;

{************************************}

      procedure Bol;
{               === }
{ This (Fortran -> Pascal) code was developed & is CopyLefted by: }
{ }
{ Han de Bruijn; Production&Services "A little bit of Physics would be  (===) }
{ TUD Computing Centre; P.O. Box 354  NO Idleness in Mathematics"(HdB) @-O^O-@ }
{ 2600 AJ  Delft; The Netherlands -- http://pubwww.tudelft.nl/~rcpshdb  #/_\# }
{ E-mail: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 278 3787. Tel: 2751.   ### }

      const NODES = 25 ;
      var xd,yd,zd,fd : row8d ;
      const zero : double = 0 ;
      var pm,term,dmin,dmax : double ;
      var k,j,i,m : byte ;

      begin
      Assign(uit,'pcbol.spl') ;
      Rewrite(uit);

      Setup_Contours;

      pm := (NODES-1)/2 ;
{ ... Coordinates: }
      for k := 1 to NODES-1 do begin
      for j := 1 to NODES-1 do begin
      for i := 1 to NODES-1 do begin
{ ... Brick values: }
      for m := 1 to 8 do begin
      xd[m] := round(i + ((m-1) mod 2))-pm ;
      yd[m] := round(j + (((m-1) div 2) mod 2))-pm ;
      zd[m] := round(k + (((m-1) div 4) mod 2))-pm ;
      term := sqr(xd[m]) + sqr(yd[m]) + sqr(zd[m]) ;
      fd[m] := sqr(10) - term ;
    { Writeln(xd[m]:5:3,' ',
              yd[m]:5:3,' ',
              zd[m]:5:3,' ',
              fd[m]:5:1,' '); }
      end;

{ ... Isosurfaces: }
      Contour_Blok(zero,xd,yd,zd,fd) ;
      end ; end ; end ;
      Close(uit);
      end;

{************************************}

procedure Kubus ;
{         ===== }

const LAYOUT = 15 ;

var w : array[1..8] of double;
    xp,yp,zp,fp : row8d ;
    i,j,k,L,m,n : byte;

function s(k : integer ; x,y,z : double) : double;
var
  a,b,c : integer;
begin
  a := (k mod 2) * 2 - 1 ;
  b := ((k div 2) mod 2) * 2 - 1 ;
  c := ((k div 4) mod 2) * 2 - 1 ;
  s := (0.5 + a*x) * (0.5 + b*y) * (0.5 + c*z) ;
end;

begin
  Assign(uit,'isocube.spl') ;
  Rewrite(uit);

  Setup_Contours;

{ ... 15 x 15 x 15 grid: }
  for i := 1 to LAYOUT-1 do
  begin
    xp[1] := 5.*(round(i-1)/round(LAYOUT-1)-0.5) ;
    xp[3] := xp[1] ; xp[5] := xp[1] ; xp[7] := xp[1] ;
    xp[2] := 5.*(round(i)/round(LAYOUT-1)-0.5) ;
    xp[4] := xp[2] ; xp[6] := xp[2] ; xp[8] := xp[2] ;

    for j := 1 to LAYOUT-1 do
    begin
      yp[1] := 5.*(round(j-1)/round(LAYOUT-1)-0.5) ;
      yp[2] := yp[1] ; yp[5] := yp[1] ; yp[6] := yp[1] ;
      yp[3] := 5.*(round(j)/round(LAYOUT-1)-0.5) ;
      yp[4] := yp[3] ; yp[7] := yp[3] ; yp[8] := yp[3] ;

      for k := 1 to LAYOUT-1 do
      begin
        zp[1] := 5.*(round(k-1)/round(LAYOUT-1)-0.5) ;
        zp[2] := zp[1] ; zp[3] := zp[1] ; zp[4] := zp[1] ;
        zp[5] := 5.*(round(k)/round(LAYOUT-1)-0.5) ;
        zp[6] := zp[5] ; zp[7] := zp[5] ; zp[8] := zp[5] ;
{ }
      { Inside / Outside function defined by: }
        for m := 1 to 8 do
        begin
          for L := 0 to 7 do
            w[L+1] := s(L,xp[m],yp[m],zp[m]) ;

          fp[m] := w[1] ;
          for L := 2 to 8 do begin
          if fp[m] > w[L] then fp[m] := w[L] ;
        end;
      end;
    {
      The function is: negative for points outside the cube,
      ................ positive for points inside the cube,
      ................ zero for points at the cube's faces. }

    { Isosurfaces: boundary / outside }
    { Contour_Blok( 0 ,xp,yp,zp,fp) ; }

      Contour_Blok( -0.65 ,xp,yp,zp,fp) ;

      end ;
    end ;
  end ;
  Close(uit);
end;

{************************************}

END.