// williams2-balls.fe
// Two Williams cells in a torus.
// Fixed version of williams.fe, with all body facets properly oriented.
// With added feature of enclosing unit balls in cells.


torus_filled
periods
1 1 0
1 -1 0
0 0  2

// for coloring contact facets according to which side they contact
define vertex attribute side integer

vertex
11 -1/4 1  1/4
12 -1/4 1 -1/4
21  1/4 1  1/4
22  1/4 1 -1/4
31  1/2 1/2  1/2
32  1/2 1/2 -1/2
51  1/2 -1/2  1/2
52  1/2 -1/2 -1/2
111 0 -1/4 1+1/4
112 0 -1/4 1-1/4
121 0  1/4 1+1/4
122 0  1/4 1-1/4

edges
10  11 12  * * *
11  11 21  * * *
12  12 22  * * *
20  21 22  * * *
21  21 31  * * *
22  22 32  * * *
31  31 11  * + *
32  32 12  * + *
41  11 51  * - *
42  12 52  * - *
51  51 21  - + *
52  52 22  - + *
110  111 112  * * *
111  111 121  * * *
112  112 122  * * *
120  121 122  * * *
121  121  32  * * +
122  122  31  * * *
131   32 111  + * -
132   31 112  + * *
141  111  52  * * +
142  112  51  * * *
151   52 121  * + -
152   51 122  * + *

faces

1  11 20 -12 -10 
2  21 31 10 -32 -22 -20 
3  41 51 20 -52 -42 -10 

11  111 120 -112 -110 
12  121 131 110 -132 -122 -120 
13  141 151 120 -152 -142 -110 

21  11 21 -122 -152 -41
22  112 122 31 41 -142 
23  11 -51 -142 -132 31
24  112 -152 51 21 132 

31  12 22 -121 -151 -42 
32  111 121 32 42 -141 
33  12 -52 -141 -131 32
34  111 -151 52 22 131 

bodies
1 -21 31 -22 32 23 -33 24 -34 -1 -2 -3 1 2 3 volume 2
2  21 -31 22 -32 -23 33 -24 34 -11 -12 -13 11 12 13 volume 2

read


// typical evolution
gogo := { 
    r; g 12; V 5; r; g 12; V 5; U;  { g 10; V} 30;
}

// draping on balls of radius 1/sqrt(2) centered at even-sum lattice points
drape := {
   local ix,iy,iz,dist,fudge,nvec;
   define nvec real[3];

   foreach vertex vv do
   {
     for ( ix := -1 ; ix <= 1 ; ix++ )
       for ( iy := -1 ; iy <= 2 ; iy++ )
         for ( iz := -1 ; iz <= 2 ; iz++ )
         { if (ix + iy + iz) imod 2 then continue;
           dist := sqrt((vv.x-ix)^2 + (vv.y-iy)^2 + (vv.z-iz)^2);
           if dist < 1/sqrt(2) then
           { fudge := 1/sqrt(2)/dist;
             vv.x := ix + fudge*(vv.x-ix);
             vv.y := iy + fudge*(vv.y-iy);
             vv.z := iz + fudge*(vv.z-iz);
             fix vv;

             // figure coloring, so can show inside sphere or outside sphere
             nvec := vv.facet[1].facet_normal;
             if nvec[1]*(vv.x-ix) + nvec[2]*(vv.y-iy) + nvec[3]*(vv.z-iz) > 0 then
               vv.side := 1
             else 
               vv.side := -1;
           }
        }
   };
   set facet color white;
   foreach facet ff where sum(ff.vertex,fixed) == 3 do
   { if sum(ff.vertex,side) > 0 then
     { ff.frontcolor := red; ff.backcolor := green; } 
     else
     { ff.frontcolor := green; ff.backcolor := red; } 
   }

}


groom := {
   unfix vertex vv where sum(vv.facet, color==white) > 0;
   g 20;
   drape;
   u; V;
}

// for good picture
run := {
  show_trans "e";  // so edges don't display; they clog things up at higher refinement
  gogo;
  groom;
  r;
  groom;
   r;
   groom 3;
   r;
   groom 3;
}

// Create picture file in postscript format.  Using view_matrix copied and pasted
// from the "print view_matrix" command after getting a good view by hand.

make_pic := {
    view_matrix := {{0.164942640162375,-0.832719935610135,-0.478958449185963,0.54929456847495},
                   {0.957591399607913,0.103858522199539,0.149204157168559,-0.746765049379614},
                   {-0.0764355846848639,-0.495803063682432,0.835683364193579,-1.22604679893956},
                   { 0, 0, 0, 1}};
    ps_colorflag on;
    ps_gridflag off;  // don't show the valence-2 edges
    full_bounding_box;  // framed in full graphics window, rather than actual bounding box
    postscript "williams-2-balls.eps";
}