// manta57adj.fe
// Adjoint of Schoen's manta genus 57 surface.

/* Commands:
   gogo - typical evolution
   showcube - show one unit cell, as on web page
   showcubelet - show 1/8 of unit cell, as on web page
   showocta - show octahedral cell, as on web page
   octa_edge - create outline edge for octahedral cell
   transforms off - show just single fundamental region
   setcolor - to color one side yellow, as in my web page.

   To turn off showing all the edges in the graphics display,
      hit the "e" key in the graphics window.
*/

parameter asize = 4.263003102112825 // shape parameter, for period killing
parameter bsize = 2.745762127998676  // shape parameter, for period killing
parameter csize = 3.492647582672936  // shape parameter, for period killing

constraint 1 // mirror plane in adjoint
formula: z - x = asize+bsize+csize+1

// Constraints for use after adjoint transformation
constraint 3 formula: y = z
constraint 4 formula: y = 0.5
constraint 5 formula: y = x
constraint 6 formula: x = 0
constraint 7 formula: x+z = 1
constraint 8 formula: z = 1

view_transform_generators 8
 -1 0 0 0   0 1 0 0    0 0 1 0    0 0 0 1  // a: x mirror
  0 1 0 0   1 0 0 0    0 0 1 0    0 0 0 1  // b: x y swap
  1 0 0 0   0 0 1 0    0 1 0 0    0 0 0 1  // c: y = z mirror
  swap_colors 0 0 -1 1  0 -1 0 1   -1 0 0 1  0 0 0 1  // d: C2 rotation
 -1 0 0 0   0 1 0 0    0 0 1 0    0 0 0 1  // e: x mirror
  1 0 0 0   0 -1 0 0   0 0 1 0    0 0 0 1  // f: y mirror
  1 0 0 0   0 1 0 0    0 0 -1 0   0 0 0 1  // g: z mirror
  1 0 0 0   0 1 0 0    0 0 -1 2   0 0 0 1  // h: z=1 mirror


vertices
1   -asize-bsize-csize-1 0 0 fixed
2   0 0 0 fixed
3   0 -1 1 fixed
4   -asize asize-1 1 fixed
5   -asize -bsize+asize-1 1+bsize fixed
6   -asize-csize  csize-bsize+asize-1  1+bsize 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 constraint 1

faces
1  1 2 3 4 5 6 

read
hessian_normal

// good evolution
gg := { refine edge where valence == 1; g 5; r; g 10; u; V;
          g5; hessian;hessian; 
          g 5; hessian; hessian;
          r; g 5; u; V; u; g 5; hessian; hessian;
          g5; hessian;hessian; 
          r; g 5; u; V; u; g 5; hessian; hessian;
          g 5; V; u; V; hessian; hessian;
          r; g 5; u; V; u; g 5; hessian; hessian;
          g 5; V; u; V; hessian; hessian;
        }

// Some distances in the adjoint
calc := {
          edge3dy := sum(edge ee where original==3,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          edge3dz := sum(edge ee where original==3,
           sum(ee.facet ff, (ff.x*ee.y-ff.y*ee.x)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge3dy - edge3dz: %g   \n",edge3dy-edge3dz;

          edge4dy := sum(edge ee where original==4,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          edge4dx := sum(edge ee where original==4,
           sum(ee.facet ff, (ff.y*ee.z-ff.z*ee.y)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge4dx - edge4dy: %g   \n",edge4dx-edge4dy;

          edge5dy := sum(edge ee where original==5,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          edge5dz := sum(edge ee where original==5,
           sum(ee.facet ff, (ff.x*ee.y-ff.y*ee.x)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge5dy - edge5dz: %g   \n",edge5dy-edge5dz;

        }


read "adjoint.cmd"

// Call this to do adjoint transformation!
adj := { unset vertex constraint 1; 
         unset edge constraint 1;
         adjoint;
       }

// Applying constraints after adjointing
frame := {
    unfix vertices; unfix edges;
    qq := vertex[1].x; set vertex x x-qq;  // vertex 1 to origin
    qq := vertex[1].y; set vertex y y-qq;
    qq := vertex[1].z; set vertex z z-qq;
    aa := maximum(max(vertex,x-y),max(vertex,y-z));
    mag := aa/.5;
    set vertex x x/mag;
    set vertex y y/mag;
    set vertex z z/mag;
    set vertex y y+.5;
    set vertex z z+1;
    foreach edge ee where original==1 do
    { set ee.vertex constraint 6; set ee constraint 6; };
    foreach edge ee where original==2 do
    { set ee constraint 3; set ee.vertex constraint 3; };
    foreach edge ee where original==3 do
    { set ee constraint 5; set ee.vertex constraint 5; };
    foreach edge ee where original==4 do
    { set ee constraint 3; set ee.vertex constraint 3; };
    foreach edge ee where original==5 do
    { set ee constraint 5; set ee.vertex constraint 5; };
    foreach edge ee where original==6 do
    { set ee constraint 7; set ee.vertex constraint 7; 
      set ee constraint 4; set ee.vertex constraint 4;
      fix ee; fix ee.vertex;
    };
}




// To get true asize after evolving after adjointing
true_asize := { printf "True asize: %20.15f\n",
   sum(edge where original == 3,length)/sum(edge where original==2,length); }
true_bsize := { printf "True bsize: %20.15f\n",
   sum(edge where original == 4,length)/sum(edge where original==2,length); }
true_csize := { printf "True csize: %20.15f\n",
   sum(edge where original == 5,length)/sum(edge where original==2,length); }
true := {true_asize;true_bsize;true_csize;}

// Set inside color to yellow
setcolor := { set facet frontcolor yellow }

// Typical evolution 
gogo := { gg; adj; frame; show_trans "R"; V; V; hessian; hessian; }

// Show 2-part manta ray
showmanta := { transform_expr "d"; show_trans "R"; }

// Show 1/8 unit cell
showcubelet := { transform_expr "cdcdcd";
                 show edge where valence == 1;  
                 show_trans "R"; 
               }

// Show unit cell
showcube := { transform_expr "gfecdcdcd"; 
              show edge where valence == 1; 
              show_trans "R"; 
            }

// Show octahedral cell
showocta    := { transform_expr "babadad"; 
                 show edge where valence <= 1;
                 show_trans "R"; 
               }
octa_edge   := { va := new_vertex(.5,.5,.5);
                  vb := new_vertex(0,0,0);
                  newe1 := new_edge(va,vb);
                  set edge[newe1] no_refine;
                  set edge[newe1] bare;
                  set edge[newe1] fixed;
                  vc := new_vertex(0,0,1);
                  newe2 := new_edge(vb,vc);
            }