// gyroid-hex.fe

// Gyroid by Bonnet rotation of hexagonal chunk of D surface.
// After loading, this file will automatically evolve the D surface,
//   and do the Bonnet rotation. 
// The user can experiment with the transform_expr command to show
// various size pieces of the gyroid.  There are 12 transform
// generators available, a-f being 4-fold rotoinversions at the
// corners of the hexagon, and g-l being C2 rotations about the 
// middles of edges.  For example:
//    transform_expr "aaa"
//    transform_expr "ghijkl"
// See the documentation on transform_expr for the full syntax.
 
// Due to lack of suitable boundary conditions, the surface is not 
// evolvable after doing the Bonnet rotation.

// replication transformations for gyroid
// Rotations and inversions about normal vectors, which are
// invariant for Bonnet rotations (except axes translated).


// Numerical coordinates along hex edge from Mathematica
#include "gyroid-bdry.inc"

// Handy array for parameterizing view transform generators
// for arbitrary Bonnet angles.
define vgtrans real[12][3]

view_transform_generators 12

// a: 4-fold rotoinversion through vertex 1
// 90 degree rotation about x axis with x inversion
swap_colors
-1 0 0 vgtrans[1][1]
 0 0 1 vgtrans[1][2] 
 0 -1 0 vgtrans[1][3]
 0 0 0 1

// b: 4-fold rotoinversion through vertex 2
// 90 degree rotation about z axis with z inversion
swap_colors
0 1 0 vgtrans[2][1]
-1 0 0 vgtrans[2][2]
 0 0 -1 vgtrans[2][3]
 0 0 0 1

// c: 4-fold rotoinversion through vertex 3
// 90 degree rotation about x axis with x inversion
swap_colors
 0 0 -1 vgtrans[3][1]
 0 -1 0 vgtrans[3][2]
 1  0 0 vgtrans[3][3]
 0 0 0 1

// d: 4-fold rotoinversion through vertex 4
// 90 degree rotation about x axis with x inversion
swap_colors
-1 0 0 vgtrans[4][1]
 0 0 1  vgtrans[4][2] 
 0 -1 0  vgtrans[4][3]
 0 0 0 1

// e: 4-fold rotoinversion through vertex 5
// 90 degree rotation about z axis with z inversion
swap_colors
0 1 0 vgtrans[5][1]
-1 0 0 vgtrans[5][2]
 0 0 -1 vgtrans[5][3]
 0 0 0 1

// f: 4-fold rotoinversion through vertex 6
// 90 degree rotation about x axis with x inversion
swap_colors
 0 0 -1 vgtrans[6][1]
 0 -1 0 vgtrans[6][2]
 1  0 0 vgtrans[6][3]
 0 0 0 1

// g: C2 rotation about midside of edge 1
0 0 1 vgtrans[7][1]
0 -1 0 vgtrans[7][2]
1 0 0 vgtrans[7][3]
0 0 0 1

// h:  C2 rotation about midside of edge 2
-1 0 0 vgtrans[8][1] 
0  0 1 vgtrans[8][2]
0  1 0 vgtrans[8][3]
0  0 0  1

// i:  C2 rotation about midside of edge 3
0  1  0  vgtrans[9][1]
1  0  0  vgtrans[9][2]
0  0 -1  vgtrans[9][3]
0  0  0  1

// j:  C2 rotation about midside of edge 4
0 0 1 vgtrans[10][1]
0 -1 0 vgtrans[10][2]
1 0 0 vgtrans[10][3]
0 0 0 1

// k:  C2 rotation about midside of edge 5
-1 0 0 vgtrans[11][1] 
0  0 1 vgtrans[11][2]
0  1 0 vgtrans[11][3]
0  0 0  1

// l:  C2 rotation about midside of edge 6
0  1  0  vgtrans[12][1]
1  0  0  vgtrans[12][2]
0  0 -1  vgtrans[12][3]
0  0  0  1

// data for D surface patch
vertices
1   1 0 0 fixed
2   1 1 0 fixed
3   0 1 0 fixed
4   0 1 1 fixed
5   0 0 1 fixed
6   1 0 1 fixed

edges
1   1 2 fixed
2   2 3 fixed
3   3 4 fixed
4   4 5 fixed
5   5 6 fixed
6   6 1 fixed

faces
1   1 2 3 4 5 6  frontcolor white backcolor yellow

read

read "adjoint.cmd"
read "adjointc.cmd"

//  Adjusting view transform generators so desired vertices are
//  fixed points.
vg_adjust := {
  for ( inx := 1 ; inx <= 6 ; inx += 1 )
  { for ( jnx := 1 ; jnx <= 3 ; jnx += 1 )
    { vgtrans[inx][jnx] := vertex[inx].__x[jnx] -
         (view_transform_generators[inx][jnx][1]*vertex[inx].__x[1] +
          view_transform_generators[inx][jnx][2]*vertex[inx].__x[2] +
          view_transform_generators[inx][jnx][3]*vertex[inx].__x[3]);
    }
  };
  for ( inx := 7 ; inx <= 12 ; inx += 1 )
  { for ( jnx := 1 ; jnx <= 3 ; jnx += 1 )
    { vgtrans[inx][jnx] := vertex[inx+1].__x[jnx] -
         (view_transform_generators[inx][jnx][1]*vertex[inx+1].__x[1] +
          view_transform_generators[inx][jnx][2]*vertex[inx+1].__x[2] +
          view_transform_generators[inx][jnx][3]*vertex[inx+1].__x[3]);
    }
  };
}

// Evolve D surface
gogoD := { r; g 5; r; g 5; hessian; r; g 5; hessian; U; g 400; vg_adjust; }

gangle := 90-atan(ellipticK(3/4)/ellipticK(1/4))*180/pi  // exact gyroid angle
            // 38.0147739891081 degrees 

// Do Bonnet rotation
bonnet := { 
  unfix vertex;
  unfix edge;
  bangle := gangle; // for default Bonnet rotation 
  adjoint;

  // center at origin
  cx := vertex[7].x; set vertex x x-cx;
  cy := vertex[7].y; set vertex y y-cy;
  cz := vertex[7].z; set vertex z z-cz;

  if valid_element(vertex[1]) then
  { // note: vertex[1] may disappear in detorus.
    // scale to unit cube 
    sc := -1/vertex[1].y;
    set vertex x x*sc;
    set vertex y y*sc;
    set vertex z z*sc;
  };

  // redefine 'g'and stuff so user doesn't accidently try to evolve
  g :::= { printf 
     "This gyroid is not evolvable due to lack of boundary constraints.\n"; 
  };
  r :::= { printf 
     "This gyroid is not evolvable due to lack of boundary constraints.\n"; 
  };
  u :::= { printf 
     "This gyroid is not evolvable due to lack of boundary constraints.\n"; 
  };
  V :::= { printf 
     "This gyroid is not evolvable due to lack of boundary constraints.\n"; 
  };
  vg_adjust;
}

detorus_sticky on; // so detorus will identify vertices with different originals

gogoD;
adjointc(0.0,7);    // sets vertex attributes so bonnet_rotate works
//bonnet;
//transform_expr "aabbcc"

//set edge color clear where valence == 2;

showq;
show_trans "R"

// Find non-redundant view transforms. Due to symmetry of the hexagon,
// there are non-identical view matrices that give the same hexagon.
// This command will sort through the current view matrices and
// print out a list of the non-redundant ones.  Output suitable
// for including as view_transforms in a datafile. Algorithm based 
// on unique copies of central vertex.

unique_transforms := {
  local center_images, unique_count, view_swaps;
  local inx,jnx,knx;

  define center_images real[transform_count][3];
  define view_list integer[transform_count];  // the unique ones 

  unique_count := 0;
  for ( inx := 1; inx <= transform_count ; inx += 1 )
  { // transform of (0,0,0) is just last column of transform matrix
    // See if close to previous images
    for ( jnx := 1 ; jnx <= unique_count ; jnx += 1 )
    { if abs(center_images[jnx][1]-view_transforms[inx][1][4]) + 
         abs(center_images[jnx][2]-view_transforms[inx][2][4]) + 
         abs(center_images[jnx][3]-view_transforms[inx][3][4]) < 0.1 then
         continue 2;
    };  
    // record new image
    unique_count += 1; 
    center_images[unique_count][1] := view_transforms[inx][1][4]; 
    center_images[unique_count][2] := view_transforms[inx][2][4]; 
    center_images[unique_count][3] := view_transforms[inx][3][4];
    view_list[unique_count] := inx;
  };
  // print out in suitable format
  printf "view_transforms %d\n", unique_count;
  for ( inx := 1 ; inx <= unique_count ; inx += 1 )
  { local ii;
    ii := view_list[inx];
    if view_transform_swap_colors[ii] then printf "swap_colors ";
    for ( jnx := 1 ; jnx <= 4; jnx += 1 )
    { for ( knx := 1 ; knx <= 4 ; knx += 1 )
        printf "%g ",view_transforms[ii][jnx][knx];
      printf "  ";
    };
    printf "\n";
  };
}

// coloring facets next to outside edges for easy transform picking
outcolor := {
  foreach edge ee where valence == 1 do set ee.facet color ee.original;
}
  
// Trying to see where symmetry is lost in bonnet rotation.
// This tests symmetry under 6-fold rotoinversion about central vertex
// with rotoinversion axis (1,1,1).
// Result: Have to be truly minimal discrete surface (not just normal-mode
// hessian) for bonnet to preserve symmetry.  But still C2 symmetry about
// outer edge midpoints not helped.

define vertex attribute mindist real
symtest := {
  midx := vertex[7].x; set vertex x x-midx;
  midy := vertex[7].y; set vertex y y-midx;
  midz := vertex[7].z; set vertex z z-midx;
  foreach vertex vv do
  { // rotate
    rx := -vv.z;
    ry := -vv.x;
    rz := -vv.y;
    // find closest vertex
    vv.mindist := 1e30;
    foreach vertex vvv do
    { dist := sqrt((vvv.x-rx)^2+(vvv.y-ry)^2+(vvv.z-rz)^2);
      if ( dist < vv.mindist ) then
         vv.mindist := dist;
    };
  } 
}

// Get hex corners exactly right, and enforce C2 symmetry about midpoint
// of each outer edge.

procedure edge_C2_fix(integer axisv,integer tailv,integer headv,real nx,real ny,real nz)
{ 
  // normalize axis vector
  mag := sqrt(nx^2 + ny^2 + nz^2);
  nx /= mag;
  ny /= mag;
  nz /= mag;

  // get a point on C2 axis 
  cx := (vertex[tailv].x + vertex[headv].x)/2;
  cy := (vertex[tailv].y + vertex[headv].y)/2;
  cz := (vertex[tailv].z + vertex[headv].z)/2;

  // get midvertex on C2 axis.
  dot_n_v := nx*(vertex[axisv].x-cx) + ny*(vertex[axisv].y-cy) + nz*(vertex[axisv].z-cz)
;
  vertex[axisv].x := cx + dot_n_v*nx;
  vertex[axisv].y := cy + dot_n_v*ny;
  vertex[axisv].z := cz + dot_n_v*nz;

  // get other vertices on edge C2 symmetric, by averaging positions
  edge_b := 0;
  foreach vertex[axisv].edge ee where valence == 1 do
  { // starting edge in each direction
    edge_a := edge_b;
    edge_b := ee.oid;
  };
  while ( edge[edge_a].vertex[2].id != headv and
          edge[edge_a].vertex[2].id != tailv ) do
  { // average positions
    va := edge[edge_a].vertex[2].id;
    vb := edge[edge_b].vertex[2].id;
    // set vb to average of vb and rotated va
    dot_n_v := nx*(vertex[va].x-cx) + ny*(vertex[va].y-cy) + nz*(vertex[va].z-cz);
    vertex[vb].x := (vertex[vb].x + vertex[va].x - 2*(vertex[va].x-cx-dot_n_v*nx))/2;    
    vertex[vb].y := (vertex[vb].y + vertex[va].y - 2*(vertex[va].y-cy-dot_n_v*ny))/2;    
    vertex[vb].z := (vertex[vb].z + vertex[va].z - 2*(vertex[va].z-cz-dot_n_v*nz))/2;    
    // now set va to the rotated value of vb
    dot_n_v := nx*(vertex[vb].x-cx) + ny*(vertex[vb].y-cy) + nz*(vertex[vb].z-cz);
    vertex[va].x := vertex[vb].x - 2*(vertex[vb].x-cx - dot_n_v*nx);
    vertex[va].y := vertex[vb].y - 2*(vertex[vb].y-cy - dot_n_v*ny);
    vertex[va].z := vertex[vb].z - 2*(vertex[vb].z-cz - dot_n_v*nz);
    // next edges
    foreach vertex[va].edge ee where valence==1 and ee.oid != -edge_a do
    { edge_a := ee.oid;
      break;
    };
    foreach vertex[vb].edge ee where valence==1 and ee.oid != -edge_b do
    { edge_b := ee.oid;
      break;
    };

  }
}

hex_tuneup := {
   // fix up corners
   for ( vinx := 1 ; vinx <= 6 ; vinx += 1 )
   { vertex[vinx].x := floor(2*vertex[vinx].x + .2)/2;
     vertex[vinx].y := floor(2*vertex[vinx].y + .2)/2;
     vertex[vinx].z := floor(2*vertex[vinx].z + .2)/2;
   };

   // Fix up sides.  Midpoint vertices are 8-13.
   edge_C2_fix(8,1,2,1,0,1);
   edge_C2_fix(9,2,3,0,1,1);
   edge_C2_fix(10,3,4,1,1,0);
   edge_C2_fix(11,4,5,1,0,1);
   edge_C2_fix(12,5,6,0,1,1);
   edge_C2_fix(13,6,1,1,1,0);

}


// Fit hexagon edge vertices to exact Mathematica coordinates.
// This assumes edge and vertex numbering as gyroid-hex.fe.
procedure exact_one_edge ( 
   integer orig_e, // original number of full hex edge
   integer xsign, integer ysign, integer zsign, // sign of coord
   integer xdim, integer ydim, integer zdim // which coord to use
)
{

   // see how many vertices along edge
   numparts := sum(edge where original==orig_e,lagrange_order);
   if numparts > max_edge_parts then
   { errprintf "exact_one_edge: Edge %d has %d parts, can only do %d.\n",
            orig_e,numparts,max_edge_parts;
     abort;
   };
   // track from tail vertex (same number as orig_e)
   tail_v := orig_e;
   this_e := 0;
   inx := 1;
   do
   {
     foreach vertex[tail_v].edge ee where 
         ee.original==orig_e and ee.oid != -this_e
     do { this_e := ee.oid; break; };

     if linear then
     { tail_v := edge[this_e].vertex[2].id;

       // fix up
     
       inx += 1;
       vertex[tail_v].x := vertex[orig_e].x 
                           + xsign*gyroid_bdrypts[numparts+1][inx][xdim];
       vertex[tail_v].y := vertex[orig_e].y 
                           + ysign*gyroid_bdrypts[numparts+1][inx][ydim];
       vertex[tail_v].z := vertex[orig_e].z 
                           + zsign*gyroid_bdrypts[numparts+1][inx][zdim];
     }
     else if quadratic then
     { tail_v := edge[this_e].vertex[2].id;
       mid_v := edge[this_e].vertex[3].id;

       // fix up
     
       inx += 1;
       vertex[mid_v].x := vertex[orig_e].x 
                           + xsign*gyroid_bdrypts[numparts+1][inx][xdim];
       vertex[mid_v].y := vertex[orig_e].y 
                           + ysign*gyroid_bdrypts[numparts+1][inx][ydim];
       vertex[mid_v].z := vertex[orig_e].z 
                           + zsign*gyroid_bdrypts[numparts+1][inx][zdim];
       inx += 1;
       vertex[tail_v].x := vertex[orig_e].x 
                           + xsign*gyroid_bdrypts[numparts+1][inx][xdim];
       vertex[tail_v].y := vertex[orig_e].y 
                           + ysign*gyroid_bdrypts[numparts+1][inx][ydim];
       vertex[tail_v].z := vertex[orig_e].z 
                           + zsign*gyroid_bdrypts[numparts+1][inx][zdim];
     }
     else // lagrange
     { 
       for ( jnx := 2 ; jnx <= lagrange_order+1 ; jnx += 1 )
       { tail_v := edge[this_e].vertex[jnx].id;

         // fix up
     
         inx += 1;
         vertex[tail_v].x := vertex[orig_e].x 
                             + xsign*gyroid_bdrypts[numparts+1][inx][xdim];
         vertex[tail_v].y := vertex[orig_e].y 
                             + ysign*gyroid_bdrypts[numparts+1][inx][ydim];
         vertex[tail_v].z := vertex[orig_e].z 
                             + zsign*gyroid_bdrypts[numparts+1][inx][zdim];
       };
     }

   } while tail_v > 7;

} // end exact_one_edge()
   

exact_edges := {
   // fix up corners
   for ( vinx := 1 ; vinx <= 6 ; vinx += 1 )
   { vertex[vinx].x := floor(2*vertex[vinx].x + .2)/2;
     vertex[vinx].y := floor(2*vertex[vinx].y + .2)/2;
     vertex[vinx].z := floor(2*vertex[vinx].z + .2)/2;
   };

   exact_one_edge(1,  1, 1,-1,  3,2,1);
   exact_one_edge(2, -1, 1,-1,  2,1,3);
   exact_one_edge(3, -1, 1, 1,  1,3,2);
   exact_one_edge(4, -1,-1, 1,  3,2,1);
   exact_one_edge(5,  1,-1, 1,  2,1,3);
   exact_one_edge(6,  1,-1,-1,  1,3,2);
   foreach edge ee where valence==1 do { fix ee; fix ee.vertex; };

   g :::= { 'g' };
   r :::= { 'r' };
   u :::= { 'u' };
   V :::= { 'V' };
}

// displaying larger piece with 12 surrounding hexagons but no overlaps.
// meant for Bonnet rotation movie.
hex13 := { transform_expr "(aa)|(bb)|(cc)|(dd)|(ee)|(ff)" }

// displaying lots and lots of hexagons.
hex_lots := { transform_expr "l(aa)|(bb)|(cc)|(dd)|(ee)|(ff)ghijk"; }

// script for making D-surface starting point for movie, 13 hexagons
// to be run immediately after loading.
dstart := { r; g 5; hessian; r; g 5; hessian; hessian; bonnet;
            exact_edges; hex13; detorus; bangle := -gangle; adjoint;
            unfix edge where valence == 2;
            unfix vertex vv where sum(vv.edge,fixed) ==  0;
            hessian; hessian;
}

text_id := -1;
label_id := -1;
screen_movie := {
   for ( frame := 0 ; frame <= 90 ; frame += 1 )
   { if frame <= 38 then ang := frame*gangle/38
     else ang := gangle + (90-gangle)/(90-38)*(frame-38); 
     bonnet_rotate(ang,7);
     angle_text := sprintf "Bonnet angle %5.2f\n",ang;
     if text_id > 0 then delete_text(text_id);
     if label_id > 0 then delete_text(label_id);
     text_id := display_text(0.1,0.1,0.03,angle_text);
     if frame == 0 then 
       label_id := display_text(0.7,0.1,0.03,"D Surface")
     else if frame == 38 then
       label_id := display_text(0.7,0.1,0.03,"Gyroid Surface")
     else if frame == 90 then 
       label_id := display_text(0.7,0.1,0.03,"P Surface");
     vg_adjust;
     recalc;
     pause;
   };
}

// Postscript frames
movie_dir := "movie/"
make_movie := {
   ps_colorflag on;
   for ( frame := 0 ; frame <= 90 ; frame += 1 )
   { if frame <= 38 then ang := frame*gangle/38
     else ang := gangle + (90-gangle)/(90-38)*(frame-38); 
     bonnet_rotate(ang,7);
     angle_text := sprintf "Bonnet angle %5.2f\n",ang;
     if text_id > 0 then delete_text(text_id);
     if label_id > 0 then delete_text(label_id);
     text_id := display_text(0.1,0.1,0.03,angle_text);
     if frame == 0 then 
       label_id := display_text(0.7,0.1,0.03,"D Surface")
     else if frame == 38 then
       label_id := display_text(0.7,0.1,0.03,"Gyroid Surface")
     else if frame == 90 then 
       label_id := display_text(0.7,0.1,0.03,"P Surface");
     vg_adjust;
     recalc;
     postscript sprintf "%sframe%03d",movie_dir,frame; 
   };
}

// Webgl movie 1: single hexagon
read "webgl-binary.cmd"
web_movie := {
   webgl_basename := "gyroid_hex";
   webgl_firstframe := 1;
   webgl_lastframe := 0;
   webgl_text_x := 10;
   webgl_text_y := 15;
   webgl_singlefile := 0;  // 1;
   webgl_indexed := 1;
   show edge where 1;
   view_matrix :=
     {{1.47646443401338,0.440336411943468,0.547136734084634,-1.23196879002074},
      {-0.501909265688248,1.55250533899715,0.104958314468604,-0.577777193888757},
      {-0.491265560042968,-0.26274132350386,1.53714902752257,-0.391571071987875},
      { 0, 0, 0, 1}};
   for ( frame := 0 ; frame <= 90 ; frame += 1 )
   { if frame <= 38 then ang := frame*gangle/38
     else ang := gangle + (90-gangle)/(90-38)*(frame-38); 
     bonnet_rotate(ang,7);
     angle_text := sprintf "Bonnet angle %5.2f\n",ang;
     if text_id > 0 then delete_text(text_id);
     if label_id > 0 then delete_text(label_id);
     text_id := display_text(0.1,0.1,0.03,angle_text);
     if frame == 0 then 
     { webgl_frame_text := sprintf"Bonnet angle %0.0f      D Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"D Surface")
     }
     else if frame == 38 then
     { webgl_frame_text := sprintf"Bonnet angle %18.15f      Gyroid Surface",gangle;
       label_id := display_text(0.7,0.1,0.03,"Gyroid Surface")
     }
     else if frame == 90 then 
     { webgl_frame_text := sprintf"Bonnet angle %3.0f      P Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"P Surface");
       webgl_lastframe := 1;
     }
     else
     {
       webgl_frame_text := sprintf"Bonnet angle %5.2f",ang;
     };
     recalc;
     if frame == 90 then webgl_lastframe := 1;
     webgl;
     webgl_firstframe := 0;
   };
}

// Webgl movie 2: three hexagons
read "webgl-binary.cmd"
web_movie_2 := {
   webgl_basename := "gyroid_hex_2";
   webgl_firstframe := 1;
   webgl_lastframe := 0;
   webgl_text_x := 10;
   webgl_text_y := 15;
   webgl_singlefile := 0;  // 1;
   webgl_indexed := 1;
   show edge where 1;
   transform_expr "aaa";
   detorus;
   adjointc(0,1);
   view_matrix :=
     {{0.727329887820938,-0.300443726303618,-0.0859782287425796,-0.727329887820937},
      {0.303750217356224,0.730861813168217,0.0156291304631574,-0.303750217356224},
      {0.0734472187966329,-0.0473500959638465,0.78678518585048,-0.0734472187966324},
      { 0, 0, 0, 1}};
   for ( frame := 0 ; frame <= 90 ; frame += 1 )
   { if frame <= 38 then ang := frame*gangle/38
     else ang := gangle + (90-gangle)/(90-38)*(frame-38); 
     bonnet_rotate(ang,1);
     angle_text := sprintf "Bonnet angle %5.2f\n",ang;
     if text_id > 0 then delete_text(text_id);
     if label_id > 0 then delete_text(label_id);
     text_id := display_text(0.1,0.1,0.03,angle_text);
     if frame == 0 then 
     { webgl_frame_text := sprintf"Bonnet angle %3.0f      D Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"D Surface")
     }
     else if frame == 38 then
     { webgl_frame_text := sprintf"Bonnet angle %18.15f      Gyroid Surface",gangle;
       label_id := display_text(0.7,0.1,0.03,"Gyroid Surface")
     }
     else if frame == 90 then 
     { webgl_frame_text := sprintf"Bonnet angle %3.0f      P Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"P Surface");
       webgl_lastframe := 1;
     }
     else
     {
       webgl_frame_text := sprintf"Bonnet angle %5.2f",ang;
     };
     recalc;
     if frame == 90 then webgl_lastframe := 1;
     webgl;
     webgl_firstframe := 0;
   };
}

// Webgl movie 3: Bigger piece.
// Getting tricky, since torus periods change with
// Bonnet angle.
read "webgl-binary.cmd"
web_movie_3 := {
   webgl_basename := "gyroid_hex_3";
   webgl_firstframe := 1;
   webgl_lastframe := 0;
   webgl_text_x := 10;
   webgl_text_y := 15;
   webgl_singlefile := 0;  // 1;
   webgl_indexed := 1;
   show edge where 1;
   transform_expr "aaa";
   detorus;
   adjointc(0,1);
   view_matrix :=
    {{0.621361369949099,0.0226937441489786,0.240498412377154,-0.442276763237615},
     {-0.00556392494008026,0.664886736203939,-0.0483644001697576,-0.30547920554705},
     {-0.241502660686568,0.0430704822473558,0.619891799336569,-0.21072981044868},
     { 0, 0, 0, 1}};
   for ( frame := 0 ; frame <= 90 ; frame += 1 )
   { replace_load datafilename;
     detorus_sticky on; // so detorus will identify vertices with different originals
     gogoD;
     adjointc(0,7);    // sets vertex attributes so bonnet_rotate works

     if frame <= 38 then ang := frame*gangle/38
     else ang := gangle + (90-gangle)/(90-38)*(frame-38); 

     bonnet_rotate(ang,7);
     vg_adjust;
     hex13;
     recalc;
     detorus; 

     angle_text := sprintf "Bonnet angle %5.2f\n",ang;
     if text_id > 0 then delete_text(text_id);
     if label_id > 0 then delete_text(label_id);
     text_id := display_text(0.1,0.1,0.03,angle_text);
     if frame == 0 then 
     { webgl_frame_text := sprintf"Bonnet angle %3.0f      D Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"D Surface")
     }
     else if frame == 38 then
     { webgl_frame_text := sprintf"Bonnet angle %18.15f      Gyroid Surface",gangle;
       label_id := display_text(0.7,0.1,0.03,"Gyroid Surface")
     }
     else if frame == 90 then 
     { webgl_frame_text := sprintf"Bonnet angle %3.0f      P Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"P Surface");
       webgl_lastframe := 1;
     }
     else
     {
       webgl_frame_text := sprintf"Bonnet angle %5.2f",ang;
     };
     recalc;
     if frame == 90 then webgl_lastframe := 1;
     webgl;
     webgl_firstframe := 0;
   };
}

// Webgl movie 4:  Really big piece. 
// Getting tricky, since torus periods change with
// Bonnet angle.
read "webgl-binary.cmd"
web_movie_4 := {
   webgl_basename := "gyroid_hex_4";
   webgl_firstframe := 1;
   webgl_lastframe := 0;
   webgl_text_x := 10;
   webgl_text_y := 15;
   webgl_singlefile := 0;  // 1;
   webgl_indexed := 1;
   show edge where 1;
   transform_expr "aaa";
   detorus;
   adjointc(0,1);
   view_matrix :=
    {{0.235344118576707,-0.0411753726435406,0.0736052616259117,-0.288361693880848},
     {0.0430823803311674,0.246259839246298,8.90999998070304e-006,-0.166221209954297},
     {-0.0725055470729925,0.0126759718386024,0.238918951491072,-0.172751390337381},
     { 0, 0, 0, 1}};
   for ( frame := 0 ; frame <= 90 ; frame += 1 )
   { replace_load datafilename;
     detorus_sticky on; // so detorus will identify vertices with different originals
     gogoD;
     adjointc(0,7);    // sets vertex attributes so bonnet_rotate works

     if frame <= 38 then ang := frame*gangle/38
     else ang := gangle + (90-gangle)/(90-38)*(frame-38); 

     bonnet_rotate(ang,7);
     vg_adjust;
     hex_lots;
     recalc;
     detorus; 

     angle_text := sprintf "Bonnet angle %5.2f\n",ang;
     if text_id > 0 then delete_text(text_id);
     if label_id > 0 then delete_text(label_id);
     text_id := display_text(0.1,0.1,0.03,angle_text);
     if frame == 0 then 
     { webgl_frame_text := sprintf"Bonnet angle %3.0f      D Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"D Surface")
     }
     else if frame == 38 then
     { webgl_frame_text := sprintf"Bonnet angle %18.15f      Gyroid Surface",gangle;
       label_id := display_text(0.7,0.1,0.03,"Gyroid Surface")
     }
     else if frame == 90 then 
     { webgl_frame_text := sprintf"Bonnet angle %3.0f      P Surface",ang;
       label_id := display_text(0.7,0.1,0.03,"P Surface");
       webgl_lastframe := 1;
     }
     else
     {
       webgl_frame_text := sprintf"Bonnet angle %5.2f",ang;
     };
     recalc;
     if frame == 90 then webgl_lastframe := 1;
     webgl;
     webgl_firstframe := 0;
   };
}