// C27Padj.fe
// Adjoint of Schoen's C27(P) surface on p. C, top.
// genus 27 surface complementary to P surface

// Programmer: Ken Brakke, brakke@susqu.edu, http://www.susqu.edu

/* Commands:
   gogo - typical evolution
   showcube - cubic unit cell, as on web page
   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.
*/

view_transform_generators 5
  0 1 0 0   1 0 0 0     0 0 1 0    0 0 0 1 // a: mirror on long curved edge
 swap_colors  // just for the next one
 -1 0 0 0   0 0 -1 0    0 -1 0 0   0 0 0 1 // b: C2 rotation on straight edge
 -1 0 0 -2  0 1 0 0     0 0 1 0    0 0 0 1 // c: mirror plane x = -1
  1 0 0 0   0 -1 0 -2   0 0 1 0    0 0 0 1 // d: mirror plane y = -1
  1 0 0 0   0 1 0 0     0 0 -1 -2  0 0 0 1 // e: mirror plane z = -1 (bottom edges)
 
// good parameters
parameter asize =  0.351182263231541  // shape parameter, for period killing
parameter bsize = 0.040487657288147  // shape parameter, for period killing

//parameter asize = 0.26  // shape parameter, for period killing
//parameter bsize = 0.23  // shape parameter, for period killing


constraint 1 // mirror plane in adjoint
formula: y = z

// Constraints for use after adjoint transformation
parameter maxy = 0
parameter minz = 0
constraint 3 formula: y = maxy
constraint 4 formula: z = minz
constraint 5 formula: y = x
constraint 6 formula: x = 0
constraint 7 formula: y+z= 0

vertices
1   0 0 0 fixed
2   -1 asize asize fixed
3   -1 asize bsize fixed
4   -1 1 bsize fixed
5   -1 1 0 fixed

edges
1  1 2 constraint 1
2  2 3 fixed
3  3 4 fixed
4  4 5 fixed
5  5 1 fixed

faces
1  -5 -4 -3 -2 -1

read
hessian_normal

// good evolution, getting lots of facets near vertex 2 cusp.
gg := { refine edge where valence == 1; g 5; r; g 10; u; V;
          refine vertex[2].edge where on_constraint 1; u; V; u; V; g5; hessian;hessian; 
          refine vertex[2].edge; u; V; u; V; g 5; hessian; hessian;
          r; g 5; u; V; u; g 5; hessian; hessian;
          refine vertex[2].edge where on_constraint 1; u; V; u; V; g5; hessian;hessian; 
          r; g 5; u; V; u; g 5; hessian; hessian;
          refine vertex[2].edge; g 5; V; u; V; u; V; u; V; hessian; hessian;
          r; g 5; u; V; u; g 5; hessian; hessian;
          refine vertex[2].edge; g 5; V; u; V; hessian; hessian;
          refine edge where original==3 or original==5; u; V;
          refine edge where original==5; u; V; g; hessian;
          hessian;
        }

// Some distances in the adjoint
calc := {
          edge1dy := sum(edge ee where original==1,
           sum(ee.facet ff, (ff.z*ee.x-ff.x*ee.z)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          edge1dz := sum(edge ee where original==1,
           sum(ee.facet ff, (ff.y*ee.x-ff.x*ee.y)/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          edge2dy := sum(edge ee where original==2,
           sum(ee.facet ff, 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/sqrt(ff.x^2+ff.y^2+ff.z^2)));
          printf " edge1dy: %g   edge1dz: %g   edge2dy: %g  edge3dz: %g \n",
              edge1dy,edge1dz,edge2dy,edge3dz;
        }

// for changing parameters
new_a := asize;
do_a  := { 
printf "do_a not working properly.\n"; return;
pp := (new_a - asize^2)/(asize-asize^2); qq := 1-pp;
           foreach vertex vv do
           { if vv on_constraint 1 then
             { unset vv constraint 1;
               vv.x := pp*x+qq*x*x;
               vv.y := pp*y+qq*y*y;
               vv.z := pp*z+qq*z*z;
               set vv constraint 1;
             } else
             { vv.x := pp*x+qq*x*x;
               vv.y := pp*y+qq*y*y;
               vv.z := pp*z+qq*z*z;
             };
           };
           asize := new_a;
         }

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;
    xoff := vertex[1].x; set vertex x x-xoff;
    yoff := vertex[1].y; set vertex y y-yoff;
    zoff := vertex[1].z; set vertex z z-zoff;
    maxy := max(vertex,y); set vertex x x/maxy; set vertex y  y/maxy; maxy := 1;
    minz := min(vertex,z); set vertex z z/(-minz); minz := -1;
    foreach edge ee where original==1 do
    { fix ee.vertex; fix ee;
      set ee.vertex constraint 6; set ee.vertex constraint 7;
      set ee constraint 6; set ee constraint 7;
    };
    foreach edge ee where original==2 do
    { set ee constraint 4; set ee.vertex constraint 4; };
    foreach edge ee where original==3 do
    { set ee constraint 3; set ee.vertex constraint 3; };
    foreach edge ee where original==4 do
    { set ee constraint 4; set ee.vertex constraint 4; };
    foreach edge ee where original==5 do
    { set ee constraint 5; set ee.vertex constraint 5; };
    
}


// for newton.cmd 
initparam := { 
   param1 := asize;
   dparam1 := 0.0001;
   target1 := 0;
}
applyparam := { new_a := param1; do_a;
                g 5; hessian; hessian; calc; 
                value1 := edge2dy;
                }

 
value2 := 0; target2 := 0; param2 := 0; dparam2 := 0; // dummies


// To get true parameters after evolving after adjointing
true_alpha := { printf "True alpha:    %20.15f\n",
   1 - sum(edge where original == 3,length)*sqrt(2)/
                sum(edge where original==5,length); }
true_beta := { printf "True beta:    %20.15f\n",
   sum(edge where original == 4,length)*sqrt(2)/
                sum(edge where original==5,length); }
true:={true_alpha;true_beta}

// for displaying full cube
showcube := { transform_expr "cdeababab"; show_trans "R";}


setcolor := { set facet backcolor yellow; }

// for nice border bands for Postscript output
read "band.cmd"
bandwidth := 0.005;
bandcolor := black;
banding := {
   edgenum := 1;
   while ( edgenum <= 5 ) do
   { set edge inband original==edgenum ; makeband;
     edgenum += 1;
   }
}


// for complete make of C27P.fe from scratch
gogo := { gg; adj;  frame; show_trans "R"; t .001; t .002; 
          hessian; hessian;
      
 }