// edge_tuber.fe
// Tubing edges.  Idea is to get the edge network you want by refining
// facets, deleting facets, deleting edges, then tubing what is left.
// Does nice junctions at vertices, and caps at valence 1 vertices.
// Same radius on all tubes.
  
/* Usage:
   Set edge_tuber_sides to desired value of sides around edge, default 6.
   Set edge_tuber_cap_zones to number of latitudinal zones wanted on edge caps,
      default 2.

   Do one of the following:
   a. To simply tube a given set of edges:
        set edge tuber_status tuber_edge_to_do where ...
        edge_tuber(radius) 
   b. To dissolve all facets and tube all edges:
        edge_tuber_naked(radius)
   c. To create dual edges and tube them:
        edge_tuber_dual(radius)

   Limitations: Only works with 3D soapfilm linear model, no torus or symmetry mode.
*/

edge_tuber_sides := 6;
edge_tuber_cap_zones := 2;
edge_tuber_maxv := edge_tuber_sides+15

// Vertices created at each vertex for each edge, stored per edge end

define edge attribute ehub integer[2][edge_tuber_maxv];  // vertex id's of hub vertices
define edge attribute ehub_count integer[2];    // how many hub vertices
define edge attribute tuber_status integer;     // edges to be tubed
tuber_edge_to_do := 1
tuber_edge_done  := 2

/* for debugging */
define facet attribute tuber_orig_edge integer  // edge the facet came from


// For adjusting tube radius after construction
define vertex attribute tuber_spot real[3]  // zero radius spot
define vertex attribute tuber_vec real[3]   // radius parametric vector

/**********************************************************************************************/

// Cap at end of edge at valence 1 vertex.
procedure edge_tuber_cap(integer v_id, real radius)
{ local sides,evec,tipv,lintelv,vx,ribe,roofe,lintele,diage,perp1,perp2;
  local inx,jnx,cangle,cc,ss,iangle,iinx,f_id;

  sides := edge_tuber_sides;

  define vx real[3];
  define evec real[3];
  define perp1 real[3];
  define perp2 real[3];
  define lintelv integer[edge_tuber_cap_zones+1][sides];
  define ribe integer[edge_tuber_cap_zones][sides];    
  define lintele integer[edge_tuber_cap_zones][sides];
  define diage integer[edge_tuber_cap_zones][sides];

  foreach vertex[v_id].edge ee where tuber_status == tuber_edge_to_do do
  {
    evec := ee.edge_vector;
    evec /= -sqrt(evec*evec);
   
    // get orthogonal basis
    if abs(evec[1]) > 0.5 then  // cross with (0,1,0)
      perp1 := { evec[3],0,-evec[1]}
    else 
      perp1 := { 0, -evec[3], evec[2]};
    perp1 /= sqrt(perp1*perp1);
    perp2 := { evec[2]*perp1[3]-evec[3]*perp1[2],
               evec[3]*perp1[1]-evec[1]*perp1[3],
               evec[1]*perp1[2]-evec[2]*perp1[1]};
    perp2 /= sqrt(perp2*perp2);
  
    // Create vertices
  
    // Tip vertex
    vx := vertex[v_id].__x + radius*evec;
    tipv := new_vertex(vx[1],vx[2],vx[3]);
    vertex[tipv].tuber_spot := vertex[v_id].__x;
    vertex[tipv].tuber_vec  := evec;
    lintelv[edge_tuber_cap_zones+1] := tipv;  // handy to have
  
    // Rings of vertices
    for ( inx := 1 ; inx <= sides ; inx++ )
    { iangle := 2*pi*inx/edge_tuber_sides;
      for ( jnx := 0 ; jnx < edge_tuber_cap_zones ; jnx++ )
      { cangle := pi/2*jnx/edge_tuber_cap_zones;
        ss := sin(cangle);
        cc := cos(cangle);
  
        vx := vertex[v_id].__x + radius*(ss*evec + cc*cos(iangle)*perp1 + cc*sin(iangle)*perp2);
        lintelv[jnx+1][inx] := new_vertex(vx[1],vx[2],vx[3]);
        vertex[lintelv[jnx+1][inx]].tuber_spot := vertex[v_id].__x;
        vertex[lintelv[jnx+1][inx]].tuber_vec := ss*evec + cc*cos(iangle)*perp1 + cc*sin(iangle)*perp2;
      }
    };
  
    // Set hub vertices
    if ee.oid > 0 then
    { ee.ehub[1] := lintelv[1];
      ee.ehub_count[1] := sides;
    }
    else
    { for ( inx := 1 ; inx <= sides ; inx++ )
        ee.ehub[2][sides-inx+1] := lintelv[1][inx];
      ee.ehub_count[2] := sides;
    };
    
    // Create edges
    for ( inx := 1 ; inx <= sides ; inx++ )
    { iinx := inx == sides ? 1 : inx+1;
      for ( jnx := 1 ; jnx <= edge_tuber_cap_zones ; jnx++ )
      {
        ribe[jnx][inx] := new_edge(lintelv[jnx][inx],lintelv[jnx+1][inx]);
        if jnx < edge_tuber_cap_zones then
          diage[jnx][inx] := new_edge(lintelv[jnx][inx],lintelv[jnx+1][iinx]);
        lintele[jnx][inx] := new_edge(lintelv[jnx][inx],lintelv[jnx][iinx]);
      }
    };
  
    // Create facets
    for ( inx := 1 ; inx <= sides ; inx++ )
    { iinx := inx == sides ? 1 : inx+1;
      for ( jnx := 1 ; jnx < edge_tuber_cap_zones ; jnx++ )
      { f_id := new_facet(lintele[jnx][inx],ribe[jnx][iinx],-diage[jnx][inx]);
        set facet[f_id] tuber_orig_edge ee.id;
        f_id := new_facet(diage[jnx][inx],-lintele[jnx+1][inx],-ribe[jnx][inx]);
        set facet[f_id] tuber_orig_edge ee.id;
      };
      f_id := new_facet(lintele[jnx][inx],ribe[jnx][iinx],-ribe[jnx][inx]);
      set facet[f_id] tuber_orig_edge ee.id;
    };
   
  } // end ee

} // end edge_tuber_cap

/**********************************************************************************************/

define vertex attribute voffset real[3] //  due to single edges

// Creating hubs at given vertex for each adjacent edge
// by finding convex set of separating planes.
// Developed from outer_vertex_tube() of dualtuber.cmd.
procedure calc_vertex_hubs(integer v_id, real radius)
{
  local plane_normals, evecs, e_oids,inx,jnx,knx,single_count,outvec;
  local mnx,ecount,pnx,hidden,found,unx,loopsize,rnx,znx,loopspot,snx; 
  local changes,sides,mminx,jeps,mray,bestj,bestdot,swangle,bestm,minx;
  local fudge, starter, maxdot, thisdot,cc,ss,mangle;
  define evecs real[20][3];  // unit edge vectors
  define plane_normals real[20][20][3]; // separating plane unit normals
  define e_oids integer[20];  // oid numbers for special edges in order found.
  
  local NOT_VISIBLE,PART_VISIBLE,WHOLE_VISIBLE;
  NOT_VISIBLE := 0;
  PART_VISIBLE := 1;
  WHOLE_VISIBLE := 2;

  // jiggle size, small enough things get identified by dtorus
  jeps := 1e-7;

  // gather edge unit vectors
  e_oids := 0;
  inx := 1; single_count := 0;
  foreach vertex[v_id].edge ee where tuber_status == tuber_edge_to_do do
  { evecs[inx] := ee.edge_vector/ee.length;
    // jiggle a bit to avoid degneracy
    evecs[inx] += { jeps*random,jeps*random,jeps*random };
    evecs[inx] /= sqrt(evecs[inx]*evecs[inx]);
    e_oids[inx] := ee.oid;
    inx++;
  };
  ecount := inx-1;

  if ecount == 0 then return;
  if ecount == 1 then { edge_tuber_cap(v_id, radius); return; };

  // calculate separating plane normals; note these point towards the inx edge, not away
  for ( inx := 1 ; inx <= ecount ; inx++ )
    for ( jnx := 1 ; jnx <= ecount ; jnx++ )
    { if inx == jnx then continue;
      plane_normals[inx][jnx] := evecs[inx]-evecs[jnx];
      plane_normals[inx][jnx] /= sqrt(plane_normals[inx][jnx]*plane_normals[inx][jnx]);
    };


  // Hub for each adjacent edge
  for ( inx := 1 ; inx <= ecount ; inx++ )
  { 
    local visible,lowray,highray,other,othercount;
    define visible integer[20];
    define lowray real[20][3];
    define highray real[20][3];
    define other integer[20];
    local invec,ray,cross,refvec,bestrefvec;
    local rays,kinx,jinx,low_neighbor,high_neighbor,raymag;
    local startjinx,other_count,mindot,bestother;
    define rays real[20][3];
    define ray real[3];
    define mray real[3];

    if ecount == 2 then 
    { othercount := 1;
    }
    else
    {
      // For each other plane, find its visibility and its neighbors
      lowray := 0; highray := 0; visible := 0; other := 0; rays := 0;

      // start with plane closest to inx edge
      bestj := 0;
      bestdot := -1;
      for ( jinx := 1 ; jinx <= ecount ; jinx++ )
      { if jinx == inx then continue;
        thisdot := evecs[jinx]*evecs[inx];
        if thisdot > bestdot then
        { bestj := jinx;
          bestdot := thisdot;
        };
      };

      kinx := bestj;
      ray := evecs[inx] + evecs[kinx];
      ray /= sqrt(ray*ray); // starter ray

      other[1] := kinx;
      othercount := 1;

      for(;;)
      {
        // now find next counterclockwise intersection ray with kinx plane
        swangle := 2*pi; 
        bestm := 0;
        for ( minx := 1 ; minx <= ecount ; minx++ )
        { if minx == inx or minx == kinx then continue;
          // intersection ray
          mray[1] := plane_normals[inx][kinx][2]*plane_normals[inx][minx][3]
                   - plane_normals[inx][kinx][3]*plane_normals[inx][minx][2];
          mray[2] := plane_normals[inx][kinx][3]*plane_normals[inx][minx][1]
                   - plane_normals[inx][kinx][1]*plane_normals[inx][minx][3];
          mray[3] := plane_normals[inx][kinx][1]*plane_normals[inx][minx][2]
                   - plane_normals[inx][kinx][2]*plane_normals[inx][minx][1];
          raymag := sqrt(mray*mray); 
          mray /= raymag;
          cc := mray*ray;
          ss := ray[1]*mray[2]*plane_normals[inx][kinx][3] +
                ray[2]*mray[3]*plane_normals[inx][kinx][1] +
                ray[3]*mray[1]*plane_normals[inx][kinx][2] -
                ray[3]*mray[2]*plane_normals[inx][kinx][1] -
                ray[1]*mray[3]*plane_normals[inx][kinx][2] -
                ray[2]*mray[1]*plane_normals[inx][kinx][3]; 
          mangle := atan2(ss,cc);
          if mangle < 0 then mangle += 2*pi;
          if mangle < swangle then
          { swangle := mangle;
            bestm := minx;
            rays[othercount] := mray;
          }
        };
        // detect loop
        if  bestm = other[1] then
          break;

        // no loop, so add to list
        othercount++;
        other[othercount] := bestm;
 
        // prepare for next go-round
        kinx := bestm;
        ray := rays[othercount-1];
      };
    }; // end ecount > 2
    
    // Create vertices and edges
    local vspot,perp1,perp2,newv,minx,cinx,fudge1,sinx,basev,mag,sidesS,zinx,zones;
    define invec real[3];
    define vspot real[3];
    define perp1 real[3];
    define perp2 real[3];
    define basev integer[edge_tuber_maxv];

    if othercount == 1 then
    { // Simple butt junction of two edges

      invec := evecs[1] + evecs[2];
      mag := sqrt(invec*invec);
      if mag < 0.00001 then
      { // straight junction
        if abs(evecs[1][1]) > 0.5 then
          invec := { -evecs[1][3], 0, evecs[1][1]}
        else
          invec := { 0, -evecs[1][3], evecs[1][2]};
        mag := sqrt(invec*invec);
      };
      perp1 := invec/mag;
      perp2 := { evecs[1][2]*perp1[3]-evecs[1][3]*perp1[2],
                 evecs[1][3]*perp1[1]-evecs[1][1]*perp1[3],
                 evecs[1][1]*perp1[2]-evecs[1][2]*perp1[1]};
      perp2 /= sqrt(perp2*perp2);
       
      sides := edge_tuber_sides;
      fudge := sqrt(2)/sqrt(1-evecs[1]*evecs[2]);
      for ( zinx := 1 ; zinx <= sides ; zinx++ )
      { vspot := vertex[v_id].__x + radius*(cos(zinx*2*pi/sides)*fudge*perp1 + sin(zinx*2*pi/sides)*perp2);
        basev[zinx] := new_vertex(vspot[1],vspot[2],vspot[3]);  
        vertex[basev[zinx]].tuber_vec := (cos(zinx*2*pi/sides)*fudge*perp1 + sin(zinx*2*pi/sides)*perp2);
        vertex[basev[zinx]].tuber_spot := vertex[v_id].__x;
      };
      foreach edge[e_oids[1]] ee do
      { 
        if ee.oid > 0 then
        { for ( zinx := 1 ; zinx <= sides ; zinx++ )
            ee.ehub[1][sides-zinx+1] := basev[zinx];
          ee.ehub_count[1] := sides;
        }
        else
        { ee.ehub[2] := basev;
          ee.ehub_count[2] := sides;
        };
      };
      foreach edge[e_oids[2]] ee do
      { 
        if ee.oid > 0 then
        { ee.ehub[1] := basev;
          ee.ehub_count[1] := sides;
        }
        else
        { for ( zinx := 1 ; zinx <= sides ; zinx++ )
            ee.ehub[2][sides-zinx+1] := basev[zinx];
          ee.ehub_count[2] := sides;
        };
      };
    }
    else 
    { // more complicated junction
      for ( minx := 1 ; minx <= othercount ; minx++ )
      { local theta1,theta2,dtheta,cc,ss,det;

        mminx := (minx == othercount) ? 1 : minx+1;
        // invec is third vector in basis of plane_normal and ray1
        cinx := other[mminx];

        invec := evecs[inx] + evecs[cinx];
        mag := sqrt(invec*invec);
        if mag < 0.00001 then
        { // straight junction
          invec := rays[minx];
          mag := sqrt(invec*invec);
        };
        invec := invec/mag;
        perp1[1] := plane_normals[inx][cinx][2]*invec[3] - plane_normals[inx][cinx][3]*invec[2];
        perp1[2] := plane_normals[inx][cinx][3]*invec[1] - plane_normals[inx][cinx][1]*invec[3];
        perp1[3] := plane_normals[inx][cinx][1]*invec[2] - plane_normals[inx][cinx][2]*invec[1];
        perp1 /= -sqrt(perp1*perp1);

        // stretch factor of intersection circle along invec
        fudge1 := sqrt(2)/sqrt(1 - (evecs[cinx]*evecs[inx]));

        // angle of ray1 
        theta1 := atan2((invec*rays[minx])/fudge1,perp1*rays[minx]);

        // angle of ray2 
        theta2 := atan2((invec*rays[mminx])/fudge1,perp1*rays[mminx]);

        dtheta := theta2 - theta1;
        if dtheta < 0 then dtheta += 2*pi;
        zones := ceil(dtheta/(2*pi)*edge_tuber_sides - 0.4978342);  // so round to nearest whole number,
                                                                    // with inexact adjustor to prevent
                                                                    // coincidences giving two different
                                                                    // zone numbers for different edges
                                                                    // at the same vertex.
        if zones == 0 then zones := 1;

        for ( jnx := 0 ; jnx < zones ; jnx++ )
        { ss := sin(theta1 + jnx/zones*dtheta);
          cc := cos(theta1 + jnx/zones*dtheta);

          vspot := vertex[v_id].__x + radius*vertex[v_id].voffset +  radius*(cc*perp1 + fudge1*ss*invec);
          newv := new_vertex(vspot[1],vspot[2],vspot[3]); 
          vertex[newv].tuber_spot := vertex[v_id].__x;
          vertex[newv].tuber_vec  := vertex[v_id].voffset + cc*perp1 + fudge1*ss*invec;
        

          foreach edge[e_oids[inx]] ee do
          { 
            if ee.oid > 0 then
            { 
              ee.ehub[1][edge_tuber_maxv - ee.ehub_count[1]] := newv;
              ee.ehub_count[1] += 1;
            }
            else
            { 
              ee.ehub_count[2] += 1;
              ee.ehub[2][ee.ehub_count[2]] := newv;
            };
          };
      
        }

      };  // end minx

      // Have to adjust hub list if tail of edge
      foreach edge[e_oids[inx]] ee do
      { 
        if ee.oid > 0 then
        {  for ( jnx := 1 ; jnx <= ee.ehub_count[1] ; jnx++ )
             ee.ehub[1][jnx] := ee.ehub[1][jnx + edge_tuber_maxv - ee.ehub_count[1]];
        }
      };
    };
 
  }; // end inx
} // end calc_vertex_hubs()

/*************************************************************************************/

// Checking for too-thick tubes, by checking for overlap of hub vertices at ends of edges.
// Returns 1 for all OK, 0 for not.
function integer edge_tuber_check()
{ local headmin,tailmax,evec,spot;
  local inx;
  define evec real[3];
 
  foreach edge ee where tuber_status == tuber_edge_to_do do
  { 
    evec := ee.edge_vector;

    tailmax := -1e30;
    for ( inx := 1; inx <= ee.ehub_count[1] ; inx++ )
    { spot := evec*vertex[ee.ehub[1][inx]].__x;
      if spot > tailmax then tailmax := spot;
    };

    headmin := 1e30;
    for ( inx := 1; inx <= ee.ehub_count[2] ; inx++ )
    { spot := evec*vertex[ee.ehub[2][inx]].__x;
      if spot < headmin then headmin := spot;
    };

    if tailmax > headmin then
    { errprintf "edge_tube_check: tube radius too small on edge %d\n",ee.id;
      return 0;
    };
  };

  return 1;

} // end edge_tuber_check()

/*************************************************************************************/

procedure make_edge_tube(integer e_id)
{ // join hubs at ends of edge, calculating explicitly the best way to join them up.
  local maxdot,bestheadhub,inx,iinx,taile,heade,newf,diagedge,evec,det;
  local tailspot,headspot,tinx,hinx,mat,sidevec,sidedot,prevedge,dotSS;
  local hinxv,tinxv;
  define heade integer[edge_tuber_maxv];
  define taile integer[edge_tuber_maxv];
  define evec real[3];
  define sidevec real[3];
  define mat real[3][3];
  define hinxv real[3];
  define tinxv real[3];


  foreach edge[e_id] ee do  // just to get handy ee name
  {
   // Check hubs have been found.
   if ee.ehub_count[1] <= 0 then
   { errprintf "edge_tuber internal error: edge %d tail hub not found.\n",e_id;
     recalc;
     abort;
   };
   if ee.ehub_count[2] <= 0 then
   { errprintf "edge_tuber internal error: edge %d head hub not found.\n",e_id;
     return;
   };

   ee.tuber_status := tuber_edge_done;

   // Find head hub vertex most aligned with first tail hub
   local evec;
   define evec real[3];
   evec := ee.edge_vector;
   maxdot := -1;
   bestheadhub := 0;
   for ( inx := 1 ; inx <= ee.ehub_count[2] ; inx++ )
   { sidevec := vertex[ee.ehub[2][inx]].__x - vertex[ee.ehub[1][1]].__x;
     sidedot := (sidevec*evec)/sqrt(sidevec*sidevec);
     if sidedot > maxdot then
     { bestheadhub := inx;
       maxdot := sidedot;
     }
   };
   if maxdot < 0.0 then
   { errprintf "edge_tuber error: can't align hubs on edge %d; edge too short for radius?\n",e_id;
     return;
   };

   // Make round-the-hub edges
   for ( inx := 1 ; inx <= ee.ehub_count[1] ; inx++ )
   { iinx := inx==ee.ehub_count[1] ? 1 : inx+1;
     taile[inx] := new_edge(ee.ehub[1][inx],ee.ehub[1][iinx]);
   };
   for ( inx := 1 ; inx <= ee.ehub_count[2] ; inx++ )
   { iinx := inx==ee.ehub_count[2] ? 1 : inx+1;
     heade[inx] := new_edge(ee.ehub[2][inx],ee.ehub[2][iinx]);
   };
     
   // Make longitudinal edges and facets, zigzagging in such a way as to keep edges
   // as close as possible to parallel to edge.  This creates an extra longitudinal edge
   // but that is eliminated in the detorus phase.
   tailspot := 1;
   headspot := bestheadhub;
   prevedge := new_edge(ee.ehub[1][tailspot],ee.ehub[2][headspot]);
   do
   {

     hinx := (headspot == ee.ehub_count[2]) ? 1 : headspot + 1;
     tinx := (tailspot == ee.ehub_count[1]) ? 1 : tailspot + 1;

     // decide on basis of which way is more parallel
     hinxv := vertex[ee.ehub[2][hinx]].__x - vertex[ee.ehub[1][tailspot]].__x;
     tinxv := vertex[ee.ehub[2][headspot]].__x - vertex[ee.ehub[1][tinx]].__x;
     if (hinxv*evec)/sqrt(hinxv*hinxv) > (tinxv*evec)/sqrt(tinxv*tinxv) then
     { // use head edge
       diagedge := new_edge(ee.ehub[1][tailspot],ee.ehub[2][hinx]);
       newf := new_facet(prevedge,heade[headspot],-diagedge);
       headspot := hinx;
       prevedge := diagedge;
     }
     else
     { // use tail edge
       diagedge := new_edge(ee.ehub[1][tinx],ee.ehub[2][headspot]);
       newf := new_facet(prevedge,-diagedge,-taile[tailspot]);
       tailspot := tinx;
       prevedge := diagedge;
     };

/* not guaranteed to make progress!

     // decide to zig or zag based on convexity of tube
     mat[1] := vertex[ee.ehub[2][hinx]].__x     - vertex[ee.ehub[1][tailspot]].__x;
     mat[2] := vertex[ee.ehub[2][headspot]].__x - vertex[ee.ehub[1][tailspot]].__x;
     mat[3] := vertex[ee.ehub[1][tinx]].__x     - vertex[ee.ehub[1][tailspot]].__x;
     det := matrix_determinant(mat);

     if det > 0 then
     { // use head edge
       diagedge := new_edge(ee.ehub[1][tailspot],ee.ehub[2][hinx]);
       newf := new_facet(prevedge,heade[headspot],-diagedge);
       headspot := hinx;
       prevedge := diagedge;
     }
     else
     { // use tail edge
       diagedge := new_edge(ee.ehub[1][tinx],ee.ehub[2][headspot]);
       newf := new_facet(prevedge,-diagedge,-taile[tailspot]);
       tailspot := tinx;
       prevedge := diagedge;
     };
*/


     set facet[newf] tuber_orig_edge e_id;

   } while ((tailspot != 1) or (headspot != bestheadhub)) ;
 } // end for

} // end make_edge_tube()

/**********************************************************************************************/

// Adjust radius of existing tubes.
procedure edge_tuber_adjust_radius(real radius)
{
  foreach vertex vv do
    vv.__x := vv.tuber_spot + radius*vv.tuber_vec;

}

/************************************************************************************/

edge_tuber_model_checks := { 

  if torus then
  { errprintf "Cannot run 'edge_tuber' command in torus mode. Do 'detorus' first.\n";
    abort;
  };

  if symmetry_group then
  { errprintf "Cannot run 'edge_tuber' command in symmetry group mode. Do 'detorus' first.\n";
    abort;
  };

  if space_dimension != 3 then
  { errprintf "The 'edge_tuber' command must be run in three-dimensional space.\n";
    abort;
  };

  if surface_dimension == 1 then
  { errprintf "The 'edge_tuber' command is not meant for the string model.\n";
    abort;
  };

  if simplex_representation then
  { errprintf "The 'edge_tuber' command is not meant for the simplex model.\n";
    abort;
  };

  if lagrange_order >= 2 then
  { errprintf "The 'edge_tuber' command is meant for the linear model, not quadratic or Lagrange.\n";
    abort;
  };

}

/*******************************************************************************/

// Basic tuber procedure.  Set edge attribute tuber_status to tuber_edge_to_do
// on edges desired to be tubed before calling. Leaves such edges with
// tuber_status set to tuber_edge_done.
procedure edge_tuber(real radius)
{ edge_tuber_model_checks;
  edge_tuber_maxv := edge_tuber_sides + 15; // room for a few extra vertices if needed
  define edge attribute ehub integer[2][edge_tuber_maxv];  // vertex id's of hub vertices
  set vertex tuber_spot __x;
  set vertex tuber_vec 0;
  set edge ehub 0;
  set edge ehub_count 0;
  foreach vertex vv do calc_vertex_hubs(vv.id,radius);
  if edge_tuber_check() == 0 then
  { dissolve vertices;
    abort;
  };

  foreach edge ee where tuber_status == tuber_edge_to_do do make_edge_tube(ee.id);
  detorus; // to identify together all the matching hub vertices and edges.
}

// Tubing all the existing edges as a naked lattice, deleting all original edges and facets.
procedure edge_tuber_naked(real radius)
{
  edge_tuber_model_checks;
  set edge tuber_status tuber_edge_to_do;
  dissolve facets;
  edge_tuber(radius);
  dissolve edges;
  dissolve vertices; 
} 

// Tubing the dual skeleton, plus special edges.
procedure edge_tuber_dual(real radius)
{
  edge_tuber_model_checks;
  set edge tuber_status tuber_edge_to_do; // mark original edges
  refine facet;  // get hub points in each facet
  edgeswap edge where tuber_status == tuber_edge_to_do and valence == 2; // dual edges
  set edge tuber_status 4 where tuber_status != tuber_edge_to_do;
  refine edge where valence != 2;
  set edge tuber_status tuber_edge_to_do where tuber_status == 0;
  set edge tuber_status 0 where tuber_status != tuber_edge_to_do;
  dissolve facets;
  edge_tuber(radius);
  dissolve edges;
  dissolve vertices; 
}


if surface_dimension != 2 then
  errprintf "\nedge_tuber: Error: tubing only works in soapfilm model.\n\n";
  
/* Usage:
   Set edge_tuber_sides to desired value of sides around edge, default 6.
   Set edge_tuber_cap_zones to number of latitudinal zones wanted on edge caps,
      default 2.

   Do one of the following:
   a. To simply tube a given set of edges:
        set edge tuber_status tuber_edge_to_do where ...
        edge_tuber(radius) 
   b. To dissolve all facets and tube all edges:
        edge_tuber_naked(radius)
   c. To create dual edges and tube them:
        edge_tuber_dual(radius)

   To adjust the tube radius afterward (no checks on overlap):
      edge_tuber_adjust_radius(radius)

   Limitations: Only works with 3D soapfilm linear model, no torus or symmetry mode.
   Will check for too-large tube radius and abort tubing.
*/