#! /usr/bin/python3
# Last edited on 2024-09-09 19:26:47 by stolfi

from math import sin, cos, tan, atan2, asin, hypot, sqrt, pi, inf, floor, sqrt
import rn, rmxn
import sys
import re
import random

import slicing_hel_povray as povray
 
prec = 4          # Decimal digits    
eps = 0.1**prec   # Coordinate quantum.

random.seed(4615)

# The program takes three parameters from the command line: 
# 
#  {pert} the coordinate perturbation amount, in units of {2*eps}.
#  {xrot} the rotation towards the {X}-axis (degrees).
#  {yrot} the rotation towards the {Y}-axis (degrees).
#
# The object is a cuboid with various protrusions and holes.

def ref_vertices(S, pock):
  # Computes the vertices of one eight of the solid. Parameters:
  
  # {S}        Side of main cube.
  
  # The value of {S} is assumed approximate. The actual dimensions
  # will be rounded to apropriate even multiples of {eps}.

  # Returns a dict with keys 'vopp', 'vpop', 'vppo', 'vppp'.
  # The values are arrays of vertices.
  # Array 'vppp' is the vertices in the positive octant.
  # Array 'vopp' has the vertices on positive quadrant of {YZ} plane.  
  # Array 'vpop' has the vertices on positive quadrant of {XZ} plane.  
  # Array 'vppo' has the vertices on positive quadrant of {XY} plane.  
  #
  # Each element is a triplet of coordinates. 
  # Does not round any coordinate.

  sys.stderr.write("ref_vertices:")
  sys.stderr.write(" S = %.*f \n" % (prec, S))
  
  # Approximate key dimensions:
  HS = S/2        # Half-side of main cube.
  Ds = HS/4       # Depth of main slot in main cube.
  Es = HS/4       # Distance in {X} from slot to edge of main cube.
  Lh = HS/2       # Length in {X} of handle.
  Wh = HS/3       # Half-width in {Y} of handle.
  Th = HS/4       # Half-height in {Z} of handle in 
  Eh = Wh/2       # Thickness of handle all around. 
  
  Xk = [ None ] * 24              
  Xk[ 0] = 1.0000 * HS      # {X} of side of main cube.
  Xk[ 1] = 0.2500 * HS      # Low {X} of slot and holes on main cube.
  Xk[ 2] = 0.7500 * HS      # High {X} of slot and holes on main cube.
  Xk[ 3] = 1.6000 * HS      # High {X} of side handle.
  Xk[ 4] = 1.2000 * HS      # Low {X} of hole in handle.
  Xk[ 5] = 1.4000 * HS      # High {X} of hole in handle. 
  Xk[ 6] = 0.4000 * HS      # Low {X} of bar in hole.
  Xk[ 7] = 0.6000 * HS      # High {X} of bar in hole.
  Xk[ 8] = 0.6250 * HS      # Some {X} of pockmarks.
  Xk[ 9] = 0.6500 * HS      # Some {X} of pockmarks.
  Xk[10] = 0.6750 * HS      # Some {X} of pockmarks.
  Xk[11] = 0.7000 * HS      # Some {X} of pockmarks.
  Xk[12] = 0.7250 * HS      # Some {X} of pockmarks.
  Xk[13] = 0.7750 * HS      # Some {X} of pockmarks.
  Xk[14] = 0.8000 * HS      # Some {X} of pockmarks.
  Xk[15] = 0.8250 * HS      # Some {X} of pockmarks.
  Xk[16] = 0.8500 * HS      # Some {X} of pockmarks.
  Xk[17] = 0.8750 * HS      # Some {X} of pockmarks.
  Xk[18] = 0.9000 * HS      # Some {X} of pockmarks.
  Xk[19] = 0.9250 * HS      # Some {X} of pockmarks.
  Xk[20] = 0.6625 * HS      # Some {X} of pockmarks.
  Xk[21] = 0.7125 * HS      # Some {X} of pockmarks.
  Xk[22] = 0.8375 * HS      # Some {X} of pockmarks.
  Xk[23] = 0.8875 * HS      # Some {X} of pockmarks.
 
  Yk = [ None ] * 4
  Yk[ 0] = 1.0000 * HS      # {Y} of side  of main cube.
  Yk[ 1] = 0.3000 * HS      # High {Y} of handle.
  Yk[ 2] = 0.1000 * HS      # High {Y} of hole in handle.
  Yk[ 3] = 0.9500 * HS      # {Y} of bottom of pockmarks.
  
  Zk = [ None ] * 22
  Zk[ 0] = 1.0000 * HS      # High {Z} of main cube.
  Zk[ 1] = 0.8000 * HS      # {Z} of bottom of slot on main cube.
  Zk[ 2] = 0.5500 * HS      # {Z} of top of hole on main cube.
  Zk[ 3] = 0.1500 * HS      # {Z} of top of handle.
  Zk[ 4] = 0.1500 * HS      # {Z} of bottom of bar in hole.
  Zk[ 5] = 0.4000 * HS      # {Z} of top of bar in hole.
  Zk[ 6] = 0.4250 * HS      # Some {Z} of pockmarks.
  Zk[ 7] = 0.4500 * HS      # Some {Z} of pockmarks.
  Zk[ 8] = 0.4750 * HS      # Some {Z} of pockmarks.
  Zk[ 9] = 0.5000 * HS      # Some {Z} of pockmarks.
  Zk[10] = 0.5250 * HS      # Some {Z} of pockmarks.
  Zk[11] = 0.5750 * HS      # Some {Z} of pockmarks.
  Zk[12] = 0.6000 * HS      # Some {Z} of pockmarks.
  Zk[13] = 0.6250 * HS      # Some {Z} of pockmarks.
  Zk[14] = 0.6500 * HS      # Some {Z} of pockmarks.
  Zk[15] = 0.6750 * HS      # Some {Z} of pockmarks.
  Zk[16] = 0.7000 * HS      # Some {Z} of pockmarks.
  Zk[17] = 0.7250 * HS      # Some {Z} of pockmarks.
  Zk[18] = 0.4625 * HS      # Some {Z} of pockmarks.
  Zk[19] = 0.5125 * HS      # Some {Z} of pockmarks.
  Zk[20] = 0.6375 * HS      # Some {Z} of pockmarks.
  Zk[21] = 0.6875 * HS      # Some {Z} of pockmarks.
  
  vopp = [ None ] *  2
  vpop = [ None ] *  0
  vppo = [ None ] *  0
  vppp = [ None ] * (15 + (28 if pock else 0))
  
  vopp[ 0] = (0,     Yk[ 0], Zk[ 0])
  vopp[ 1] = (0,     Yk[ 0], Zk[ 2])

  vppp[ 0] = (Xk[ 1], Yk[ 0], Zk[ 0])
  vppp[ 1] = (Xk[ 1], Yk[ 0], Zk[ 1])
  vppp[ 2] = (Xk[ 2], Yk[ 0], Zk[ 1])
  vppp[ 3] = (Xk[ 2], Yk[ 0], Zk[ 0])
  vppp[ 4] = (Xk[ 0], Yk[ 0], Zk[ 0])
  vppp[ 5] = (Xk[ 1], Yk[ 0], Zk[ 2])
  vppp[ 6] = (Xk[ 2], Yk[ 0], Zk[ 2])
  vppp[ 7] = (Xk[ 0], Yk[ 1], Zk[ 3])
  vppp[ 8] = (Xk[ 3], Yk[ 1], Zk[ 3])
  vppp[ 9] = (Xk[ 4], Yk[ 2], Zk[ 3])
  vppp[10] = (Xk[ 5], Yk[ 2], Zk[ 3])

  #  Bars in holes
  vppp[11] = (Xk[ 6], Yk[ 0], Zk[ 5])
  vppp[12] = (Xk[ 7], Yk[ 0], Zk[ 5])
  vppp[13] = (Xk[ 6], Yk[ 0], Zk[ 4])
  vppp[14] = (Xk[ 7], Yk[ 0], Zk[ 4])
  
  if pock:
    # Pockmark top corners:
    vppp[15] = (Xk[ 2], Yk[ 0], Zk[12])
    vppp[16] = (Xk[11], Yk[ 0], Zk[12])
    vppp[17] = (Xk[11], Yk[ 0], Zk[14])
    vppp[18] = (Xk[10], Yk[ 0], Zk[12])
    vppp[19] = (Xk[ 8], Yk[ 0], Zk[13])
    vppp[20] = (Xk[ 9], Yk[ 0], Zk[15])
    vppp[21] = (Xk[10], Yk[ 0], Zk[16])
    vppp[22] = (Xk[12], Yk[ 0], Zk[17])
    vppp[23] = (Xk[ 2], Yk[ 0], Zk[15])
    vppp[24] = (Xk[ 2], Yk[ 0], Zk[14])
    vppp[25] = (Xk[14], Yk[ 0], Zk[12])
    vppp[26] = (Xk[14], Yk[ 0], Zk[ 2])
    vppp[27] = (Xk[16], Yk[ 0], Zk[ 2])
    vppp[28] = (Xk[16], Yk[ 0], Zk[ 9])
    vppp[29] = (Xk[17], Yk[ 0], Zk[ 2])
    vppp[30] = (Xk[19], Yk[ 0], Zk[10])
    vppp[31] = (Xk[18], Yk[ 0], Zk[ 8])
    vppp[32] = (Xk[17], Yk[ 0], Zk[ 7])
    vppp[33] = (Xk[15], Yk[ 0], Zk[ 6])
    vppp[34] = (Xk[14], Yk[ 0], Zk[ 8])
    vppp[35] = (Xk[14], Yk[ 0], Zk[ 9])

    # Pockmark bottoms:
    vppp[36] = (Xk[15], Yk[ 3], Zk[10])
    vppp[37] = (Xk[12], Yk[ 3], Zk[13])
    vppp[38] = (Xk[22], Yk[ 3], Zk[18])
    vppp[39] = (Xk[23], Yk[ 3], Zk[19])
    vppp[40] = (Xk[20], Yk[ 3], Zk[20])
    vppp[41] = (Xk[21], Yk[ 3], Zk[21])
    vppp[42] = (Xk[13], Yk[ 3], Zk[11])

  return { 'opp': vopp, 'pop': vpop, 'ppo': vppo, 'ppp': vppp }
 
def all_vertices(vref):
  # The dict {vref} must be a {dict} that gives the vertices 
  # on the positive octant (key 'ppp') and on the planes 
  # {YZ} ('opp'), {XZ} ('pop'), and {XY} ('ppo').
  #
  # Returns a dict {vtot} completed with vertices 
  # on other octants and quadrants of those planes,
  # with keys 'mpp', 'mpm', 'omp', etc.  Each vertex is a quadruple 
  # where the first three elements are the coordinates, 
  # followed by the OBJ index of the vertex (from 1).
  #
  # Also returns the number of vertices created {Nv}.

  vtot = {}
  
  Nv = 0 # Number of vertices generated so far.
  
  for sx in +1, 0, -1:
    for sy in +1, 0, -1:
      for sz in +1, 0, -1:
        key_src = f"{'pop'[sx+1]}{'pop'[sy+1]}{'pop'[sz+1]}"
        v_src = vref.get(key_src, None)
        if v_src != None and len(v_src) != 0:
          Nv_src = len(v_src)
          key_dst = f"{'mop'[sx+1]}{'mop'[sy+1]}{'mop'[sz+1]}"
          v_dst = [ None ] * Nv_src
          for iv in range(Nv_src):
            Nv += 1
            vi_src = v_src[iv]
            vi_dst = ( vi_src[0]*sx, vi_src[1]*sy, vi_src[2]*sz, Nv, ) 
            v_dst[iv] = vi_dst
          vtot[key_dst] = v_dst
  return vtot, Nv
  
def perturb_and_round_vertices(vtot, rmat,pert):
  # Multiplies eack coord vector of each vertex of {vtot} (as row) by the 3x3 matrix {rmat}.
  # Then adds a random amount in {[-pert _ pert]} to each coordinate.
  # Then rounds each coordinate to an even multiple of {eps}.
  
  for key,v in vtot.items():
    for iv in range(len(v)):
      vi = v[iv]
      ci = rmxn.map_row(vi[:3], rmat)
      v[iv] = tuple([ pertround(ci[j],pert) for j in range(3) ] + [ vi[3], ])
  return None

def pertround(c,pert):
  # Adds a random amount in {[-pert _ pert]*2*eps} to {c}, then 
  # rounds to an even multiple of {eps}.
  
  if pert != 0: 
    d = random.uniform(-pert*2*eps, pert*2*eps)
    c = 2*eps * floor(c/(2*eps) + 0.5)
  return c
  
def write_vertices_OBJ(wro, vtot):
  # Writes the vertices of the object to the file {wro} in OBJ format.
  # 
  # Expect {vtot} to be a dict with keys 'mmm', 'mmo', 'mmp', 'mom', ..., 'ppm',
  # 'ppo,' 'ppp' with the vertices in each 'octant' of {\RR^3},
  # Writes the vertices to {wro} in OBJ format.
  # 
  # Assumes that each vertex in {vtot} has three {float}
  # coordinates and a fourth integer element that is the OBJ index.
  # Assumes that each coordinate is rounded to an even multiple of
  # {eps}.
    
  debug_vertices = True
  
  Nv = 0
  
  if debug_vertices:
    sys.stderr.write("writing vertices...\n")

  def Prtv(p, lab):
    # Writes {p[0..2]} to {wro} as the coordinates of a new vertex.
    # Expects that the fourth element of {p} is the OBJ index of the vertex
    # and vertices are printed in the order of that index.
    # The {lab} is the vertex label "{key}:{index}" in {vtot}.
    nonlocal Nv
    Nv += 1;
    assert Nv == p[3]
    # Writes {p} to {wro}: 
    wro.write("v");
    for i in range(3): wro.write(" %.*f" % (prec, p[i]))
    wro.write(" # %s\n" % lab);
    if debug_vertices:
      sys.stderr.write("  V[%d] = %s = (" % (p[3], lab))
      for i in range(3): sys.stderr.write(" %.*f" % (prec, p[i]))
      sys.stderr.write(" )\n")

  for key,v in vtot.items():
    for iv in range(len(v)):
      lab = f"{key}:{iv}"
      Prtv(v[iv], lab)

  wro.write("\n")

  return None
  
def write_faces_OBJ(wro, vtot, pock):
  # Writes the faces of the object to OBJ file {wro}, assuming that
  # the vertices are in the dict {vtot}.
  
  # Also returns tables {Fnrm} of normals to the faces and {Fbar} of their
  # barycenters.
  #
  # Also returns tables {E,G} of edges:
  #
  #   {E} Non-ghost edges of the mesh.
  #   {G} Ghost edges of the mesh.
  #
  # Each element of these arrays is a pair {(iorg,idst)} of OBJ indices
  # of the endpoint vertices. Each array {Fnrm,Fbar,E,G} is indexed
  # {[1..N]}; element 0 is not used ({None}).
  
  # Expected face, and edge counts perpendicular or parallel to each axis, and
  # general oblique:
  
  #        main shape      
  #        -----------------
  #           X     Y     Z
  #        -----------------       
  Nf_exp =   28 +  18 +  26   # Expected face counts (perp).
  Ne_exp =   64 +  60 +  60   # Expected non-ghost edge count (para).
  Ng_exp =    2 +   6 +   4   # Expected ghost edge count (perp).
  
  if pock:
    #           ----------------------
    #              X     Y     Z  oblq
    #           ----------------------
    Nf_exp +=      0 +   0 +   0 + 224
    Ne_exp +=     48 +   0 +  48 + 288
    Ng_exp +=      0 +   0 +   0 +   0
    
  # Variables that are global to the nested functions:

  E = [ None ]
  G = [ None ]
  Fbar = [ None ]
  Fnrm = [ None ]

  debug_face = False
  Fp = []    # Vertices of current face, indexed from 0.
  
  def Sved(p_org, p_dst, ghost):
    # Saves the edge from {p_org} to {p_dst} in {E} or {G} depending on
    # {ghost}, if traversed in increasing vertex index sense. 

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    kv_org = p_org[3]
    kv_dst = p_dst[3]
    # Avoid duplicate edges (even ghost ones):
    if kv_org < kv_dst:
      if ghost:
        G.append((kv_org, kv_dst))
      else:
        E.append((kv_org, kv_dst))
    return None

  def Ptit(tit):
    # Title for a face or set of faces.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    wro.write("\n");
    wro.write("# %s\n" % tit);
    sys.stderr.write("%s \n" % (tit));
    return None
  
  def Bof(debug = False):
    #  Start of a new face.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    debug_face = debug
    assert len(Fp) == 0
    return None
  
  def Prtv(key, iv, ghost):
    # Adds vertex {vtot[key][iv]} to the current face, augmented with 
    # its "{key}:{iv}" label.  
    #
    # If {ghost} is {True}, assumes that the edge from previous
    # vertex is a ghost one.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    p = vtot[key][iv]
    assert len(p) == 4
    # Saves vertex in list {Fp} of face vertices, with label:
    lab = f"{key}:{iv:02}"
    Fp.append(p + ( lab, ));
    if len(Fp) >= 2: Sved(Fp[-2], Fp[-1], ghost)
    return None
    
  def Eof(ghost, order):
    # End of face. If {ghost}, assumes that 
    # closing edge is a ghost one.
    # Prints in given or reverse order depending on {order} is {+1} or {-1}.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face

    if debug_face:
      sys.stderr.write("  writing face %d (%d corners)\n" % (len(Fnrm), len(Fp))) 
    assert len(Fp) >= 3
    Sved(Fp[-1], Fp[0], ghost)

    # Write the face:
    wro.write("f")
    deg = len(Fp)
    for rv in range(deg):
      jv = rv if order == +1 else deg-1-rv
      p = Fp[jv]
      assert len(p) == 5
      kv = p[3]
      lab = p[4]
      if debug_face:
        sys.stderr.write("    v%04d = %s = ( %9.*f %9.*f %9.*f )\n" % (kv, lab, prec, p[0], prec, p[1], prec, p[2]));
      wro.write(" %d" % kv);
    wro.write("\n");
    
    # Compute normal, barycenter:
    nrm, bar = face_normal_and_barycenter(Fp)
    Fnrm.append(nrm)
    Fbar.append(bar)
    Fp = []; 
    return None
    
  # Here are the faces:
  
  debug_faces = True

  Ptit(f"faces with normal ±X")
  for sx in +1, -1:
    tx = 'mop'[sx+1]

    Ptit(f"  faces in halfspace sx = {sx:+d}")
    
    # Vert face of cube:
    Bof(debug_faces)
    Prtv(tx+'pp',  4, False)
    
    Prtv(tx+'mp',  4, False)
    Prtv(tx+'mm',  4, False)
    Prtv(tx+'pm',  4, False)
    Prtv(tx+'pp',  4, False)
    
    Prtv(tx+'pp',  7, True)
    Prtv(tx+'pm',  7, False)
    Prtv(tx+'mm',  7, False)
    Prtv(tx+'mp',  7, False)
    Prtv(tx+'pp',  7, False)
    Eof(True, sx)
    
    # Centrad vert wall of hole on cube:
    Bof(debug_faces)
    Prtv(tx+'pp',  5, False)
    Prtv(tx+'mp',  5, False)
    Prtv(tx+'mm',  5, False)
    Prtv(tx+'pm',  5, False)
    Eof(False, sx)
    
    # Distad vert wall of hole on cube:
    Bof(debug_faces)
    Prtv(tx+'pp',  6, False)
    Prtv(tx+'pm',  6, False)
    Prtv(tx+'mm',  6, False)
    Prtv(tx+'mp',  6, False)
    Eof(False, sx)
    
    # Centrad vert wall of hole on handle:
    Bof(debug_faces)
    Prtv(tx+'pp',  9, False)
    Prtv(tx+'mp',  9, False)
    Prtv(tx+'mm',  9, False)
    Prtv(tx+'pm',  9, False)
    Eof(False, sx)
    
    # Distad vert wall of hole on handle:
    Bof(debug_faces)
    Prtv(tx+'pp', 10, False)
    Prtv(tx+'pm', 10, False)
    Prtv(tx+'mm', 10, False)
    Prtv(tx+'mp', 10, False)
    Eof(False, sx)
    
    # Extreme vert wall of handle:
    Bof(debug_faces)
    Prtv(tx+'pp',  8, False)
    Prtv(tx+'mp',  8, False) 
    Prtv(tx+'mm',  8, False)
    Prtv(tx+'pm',  8, False)
    Eof(False, sx)

    for sz in +1, -1:
      tz = 'mop'[sz+1]

      Ptit(f"    faces in quadrant sx = {sx:+d}, sz = {sz:+d}")
      
      # Centrad vert wall of slot on face of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  0, False)
      Prtv(tx+'m'+tz,  0, False)
      Prtv(tx+'m'+tz,  1, False)
      Prtv(tx+'p'+tz,  1, False)
      Eof(False, sx*sz)
      
      # Distad vert wall of slot on {±Z} face of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  2, False)
      Prtv(tx+'m'+tz,  2, False)
      Prtv(tx+'m'+tz,  3, False)
      Prtv(tx+'p'+tz,  3, False)
      Eof(False, sx*sz)
       
      # Distad vert wall of bar in hole of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz, 12, False)
      Prtv(tx+'m'+tz, 12, False)
      Prtv(tx+'m'+tz, 14, False)
      Prtv(tx+'p'+tz, 14, False)
      Eof(False, sx*sz)
      
      # Centrad vert wall of bar in hole of cube:
      Bof(debug_faces)
      Prtv(tx+'m'+tz, 11, False)
      Prtv(tx+'p'+tz, 11, False)
      Prtv(tx+'p'+tz, 13, False)
      Prtv(tx+'m'+tz, 13, False)
      Eof(False, sx*sz)

  Ptit(f"faces with normal ±Y")

  for sy in +1, -1:
    ty = 'mop'[sy+1]
    
    Ptit(f"  faces in halfpspace sy = {sy:+d}")
    
    # Vert face of main cube:
    Ptit(f"    big complicated face with holes and pinch corners:")
    
    Prtv('o'+ty+'m',  0, False)
    Prtv('o'+ty+'m',  1, True)
    
    Prtv('p'+ty+'m',  5, True)
    Prtv('p'+ty+'m',  6, False)

    if pock:
      Prtv('p'+ty+'m', 15, False)
      Prtv('p'+ty+'m', 16, False)
      Prtv('p'+ty+'m', 17, False)
      Prtv('p'+ty+'m', 18, False)
      Prtv('p'+ty+'m', 19, False)
      Prtv('p'+ty+'m', 20, False)
      Prtv('p'+ty+'m', 17, False)
      Prtv('p'+ty+'m', 21, False)
      Prtv('p'+ty+'m', 22, False)
      Prtv('p'+ty+'m', 23, False)
      Prtv('p'+ty+'m', 17, False)
      Prtv('p'+ty+'m', 24, False)
      Prtv('p'+ty+'m', 15, False)
      Prtv('p'+ty+'m', 25, False)
      Prtv('p'+ty+'m', 26, False)
      Prtv('p'+ty+'m', 27, False)
      Prtv('p'+ty+'m', 28, False)
      Prtv('p'+ty+'m', 29, False)
      Prtv('p'+ty+'m', 30, False)
      Prtv('p'+ty+'m', 31, False)
      Prtv('p'+ty+'m', 28, False)
      Prtv('p'+ty+'m', 32, False)
      Prtv('p'+ty+'m', 33, False)
      Prtv('p'+ty+'m', 34, False)
      Prtv('p'+ty+'m', 28, False)
      Prtv('p'+ty+'m', 35, False)
      Prtv('p'+ty+'m', 26, False)
      Prtv('p'+ty+'m',  6, False)

    Prtv('p'+ty+'p',  6, False)

    if pock:
      Prtv('p'+ty+'p', 26, False)
      Prtv('p'+ty+'p', 35, False)
      Prtv('p'+ty+'p', 28, False)
      Prtv('p'+ty+'p', 34, False)
      Prtv('p'+ty+'p', 33, False)
      Prtv('p'+ty+'p', 32, False)
      Prtv('p'+ty+'p', 28, False)
      Prtv('p'+ty+'p', 31, False)
      Prtv('p'+ty+'p', 30, False)
      Prtv('p'+ty+'p', 29, False)
      Prtv('p'+ty+'p', 28, False)
      Prtv('p'+ty+'p', 27, False)
      Prtv('p'+ty+'p', 26, False)
      Prtv('p'+ty+'p', 25, False)
      Prtv('p'+ty+'p', 15, False)
      Prtv('p'+ty+'p', 24, False)
      Prtv('p'+ty+'p', 17, False)
      Prtv('p'+ty+'p', 23, False)
      Prtv('p'+ty+'p', 22, False)
      Prtv('p'+ty+'p', 21, False)
      Prtv('p'+ty+'p', 17, False)
      Prtv('p'+ty+'p', 20, False)
      Prtv('p'+ty+'p', 19, False)
      Prtv('p'+ty+'p', 18, False)
      Prtv('p'+ty+'p', 17, False)
      Prtv('p'+ty+'p', 16, False)
      Prtv('p'+ty+'p', 15, False)
      Prtv('p'+ty+'p',  6, False)

    Prtv('p'+ty+'p',  5, False)
    Prtv('p'+ty+'m',  5, False)
    
    Prtv('o'+ty+'m',  1, True)
    
    Prtv('m'+ty+'m',  5, True)
    Prtv('m'+ty+'p',  5, False)
    
    Prtv('m'+ty+'p',  6, False)

    if pock:
      Prtv('m'+ty+'p', 15, False)
      Prtv('m'+ty+'p', 16, False)
      Prtv('m'+ty+'p', 17, False)
      Prtv('m'+ty+'p', 18, False)
      Prtv('m'+ty+'p', 19, False)
      Prtv('m'+ty+'p', 20, False)
      Prtv('m'+ty+'p', 17, False)
      Prtv('m'+ty+'p', 21, False)
      Prtv('m'+ty+'p', 22, False)
      Prtv('m'+ty+'p', 23, False)
      Prtv('m'+ty+'p', 17, False)
      Prtv('m'+ty+'p', 24, False)
      Prtv('m'+ty+'p', 15, False)
      Prtv('m'+ty+'p', 25, False)
      Prtv('m'+ty+'p', 26, False)
      Prtv('m'+ty+'p', 27, False)
      Prtv('m'+ty+'p', 28, False)
      Prtv('m'+ty+'p', 29, False)
      Prtv('m'+ty+'p', 30, False)
      Prtv('m'+ty+'p', 31, False)
      Prtv('m'+ty+'p', 28, False)
      Prtv('m'+ty+'p', 32, False)
      Prtv('m'+ty+'p', 33, False)
      Prtv('m'+ty+'p', 34, False)
      Prtv('m'+ty+'p', 28, False)
      Prtv('m'+ty+'p', 35, False)
      Prtv('m'+ty+'p', 26, False)
      Prtv('m'+ty+'p',  6, False)

    Prtv('m'+ty+'m',  6, False)

    if pock:
      Prtv('m'+ty+'m', 26, False)
      Prtv('m'+ty+'m', 35, False)
      Prtv('m'+ty+'m', 28, False)
      Prtv('m'+ty+'m', 34, False)
      Prtv('m'+ty+'m', 33, False)
      Prtv('m'+ty+'m', 32, False)
      Prtv('m'+ty+'m', 28, False)
      Prtv('m'+ty+'m', 31, False)
      Prtv('m'+ty+'m', 30, False)
      Prtv('m'+ty+'m', 29, False)
      Prtv('m'+ty+'m', 28, False)
      Prtv('m'+ty+'m', 27, False)
      Prtv('m'+ty+'m', 26, False)
      Prtv('m'+ty+'m', 25, False)
      Prtv('m'+ty+'m', 15, False)
      Prtv('m'+ty+'m', 24, False)
      Prtv('m'+ty+'m', 17, False)
      Prtv('m'+ty+'m', 23, False)
      Prtv('m'+ty+'m', 22, False)
      Prtv('m'+ty+'m', 21, False)
      Prtv('m'+ty+'m', 17, False)
      Prtv('m'+ty+'m', 20, False)
      Prtv('m'+ty+'m', 19, False)
      Prtv('m'+ty+'m', 18, False)
      Prtv('m'+ty+'m', 17, False)
      Prtv('m'+ty+'m', 16, False)
      Prtv('m'+ty+'m', 15, False)
      Prtv('m'+ty+'m',  6, False)

    Prtv('m'+ty+'m',  5, False)
    
    Prtv('o'+ty+'m',  1, True)
    Prtv('o'+ty+'m',  0, True)
    
    Prtv('m'+ty+'m',  0, False)
    Prtv('m'+ty+'m',  1, False)
    Prtv('m'+ty+'m',  2, False)
    Prtv('m'+ty+'m',  3, False)
    Prtv('m'+ty+'m',  4, False)

    Prtv('m'+ty+'p',  4, False)
    Prtv('m'+ty+'p',  3, False)
    Prtv('m'+ty+'p',  2, False)
    Prtv('m'+ty+'p',  1, False)
    Prtv('m'+ty+'p',  0, False)
    
    Prtv('o'+ty+'p',  0, False)
    
    Prtv('p'+ty+'p',  0, False)
    Prtv('p'+ty+'p',  1, False)
    Prtv('p'+ty+'p',  2, False)
    Prtv('p'+ty+'p',  3, False)
    Prtv('p'+ty+'p',  4, False)

    Prtv('p'+ty+'m',  4, False)
    Prtv('p'+ty+'m',  3, False)
    Prtv('p'+ty+'m',  2, False)
    Prtv('p'+ty+'m',  1, False)
    Prtv('p'+ty+'m',  0, False)

    Eof(False, sy)
     
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sx = {sx:+d} sy = {sy:+d}")
    
      # Vert outer face of handle:
      Bof(debug_faces)
      Prtv(tx+ty+'p',  7, False)
      Prtv(tx+ty+'p',  8, False)
      Prtv(tx+ty+'m',  8, False)
      Prtv(tx+ty+'m',  7, False)
      Eof(False, sx*sy)
    
      # Vert inner face of handle:
      Bof(debug_faces)
      Prtv(tx+ty+'p', 10, False)
      Prtv(tx+ty+'p',  9, False)
      Prtv(tx+ty+'m',  9, False)
      Prtv(tx+ty+'m', 10, False)
      Eof(False, sx*sy)
      
      for sz in +1, -1:
        tz = 'mop'[sz+1]

        Ptit(f"      faces in octant sx = {sx:+d} sy = {sy:+d} sz = {sz:+d}")

        # Tip of bar in hole:
        Bof(debug_faces)
        Prtv(tx+ty+tz, 11, False)
        Prtv(tx+ty+tz, 12, False)
        Prtv(tx+ty+tz, 14, False)
        Prtv(tx+ty+tz, 13, False)
        Eof(False, sx*sy*sz)

        if pock:
          Ptit(f"      pockmark faces in octant sx = {sx:+d} sy = {sy:+d} sz = {sz:+d}")
          # Each tuple below is center vertex and perimeter vertices in CCW order.
          Pockmarks = (
            ( 36, 26, 35, 28, 27, ),
            ( 37, 15, 24, 17, 16, ),
            ( 38, 28, 34, 33, 32, ),
            ( 39, 28, 31, 30, 29, ),
            ( 40, 17, 20, 19, 18, ),
            ( 41, 17, 23, 22, 21, ),
            ( 42,  6, 26, 25, 15, ),
          )
          for Pck in Pockmarks:
            ivt = Pck[0]
            for ka in range(4):
              kb = (ka + 1) % 4
              iva = Pck[ka+1]
              ivb = Pck[kb+1]
              Bof(debug_faces)
              Prtv(tx+ty+tz, ivt, False)
              Prtv(tx+ty+tz, iva, False)
              Prtv(tx+ty+tz, ivb, False)
              Eof(False, sx*sy*sz)

  Ptit(f"faces with normal ±Z")

  for sz in +1, -1:
    tz = 'mop'[sz+1]
    
    Ptit(f"  faces in halfspace sz = {sz:+d}")
    
    # Centrad horz face of main cube:
    Bof(debug_faces)
    Prtv('op'+tz,  0, False)
    Prtv('mp'+tz,  0, False)
    Prtv('mm'+tz,  0, False)
    Prtv('om'+tz,  0, False)
    Prtv('pm'+tz,  0, False)
    Prtv('pp'+tz,  0, False)
    Eof(False, sz)
    
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sz = {sz:+d}, sx = {sx:+d}")

      # Distad horz face of main cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  3, False)
      Prtv(tx+'m'+tz,  3, False)
      Prtv(tx+'m'+tz,  4, False)
      Prtv(tx+'p'+tz,  4, False)
      Eof(False, sx*sz)
     
      # Horz floor of slot on top of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  1, False)
      Prtv(tx+'m'+tz,  1, False)
      Prtv(tx+'m'+tz,  2, False)
      Prtv(tx+'p'+tz,  2, False)
      Eof(False, sx*sz)
      
      # Horz ceil of hole on main cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  5, False)
      Prtv(tx+'p'+tz,  6, False)
      Prtv(tx+'m'+tz,  6, False)
      Prtv(tx+'m'+tz,  5, False)
      Eof(False, sx*sz)
      
      # Top horz face of bar in hole:
      Bof(debug_faces)
      Prtv(tx+'p'+tz, 12, False)
      Prtv(tx+'p'+tz, 11, False)
      Prtv(tx+'m'+tz, 11, False)
      Prtv(tx+'m'+tz, 12, False)
      Eof(False, sx*sz)
      
      # Bot horz face of bar in hole:
      Bof(debug_faces)
      Prtv(tx+'p'+tz, 13, False)
      Prtv(tx+'p'+tz, 14, False)
      Prtv(tx+'m'+tz, 14, False)
      Prtv(tx+'m'+tz, 13, False)
      Eof(False, sx*sz)
        
      # Horz face on handle:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  7, False)
      Prtv(tx+'m'+tz,  7, False)
      Prtv(tx+'m'+tz,  8, False)
      Prtv(tx+'p'+tz,  8, False)
      Prtv(tx+'p'+tz,  7, False)
      Prtv(tx+'p'+tz,  9, True)
      Prtv(tx+'p'+tz, 10, False)
      Prtv(tx+'m'+tz, 10, False)
      Prtv(tx+'m'+tz,  9, False)
      Prtv(tx+'p'+tz,  9, False)
      Eof(True, sx*sz)
  
  sys.stderr.write("wrote %d faces (expected %d)\n" % (len(Fnrm)-1, Nf_exp))
  assert len(Fnrm) == Nf_exp + 1
  assert len(Fbar) == Nf_exp + 1

  sys.stderr.write("found %d non-ghost edges (expected %d)\n" % (len(E)-1, Ne_exp))
  assert len(E) == Ne_exp + 1

  sys.stderr.write("found %d ghost edges (expected %d)\n" % (len(G)-1, Ng_exp))
  assert len(G) == Ng_exp + 1

  return Fnrm, Fbar, E, G
  
def face_normal_and_barycenter(Fp):
  # Compute normal {nrm} and barycenter {bar}
  # given a list of vertices {Fp[0..deg-1]}
  deg = len(Fp)
  assert deg >= 3
  po = Fp[0][:3]
  pa = Fp[1][:3]
  pb = Fp[2][:3]
  u = rn.sub(pa,po)
  v = rn.sub(pb,po)
  nrm, sz = rn.dir(rn.cross3d(u, v))

  bar = (0, 0, 0)
  for jv in range(deg):
    pj = Fp[jv][:3]
    bar = rn.mix(1.0, bar, 1.0/deg, pj)
    
  # Check planarity:
  for jv in range(deg):
    pj = Fp[jv][:3]
    u = rn.sub(pj, bar)
    s = rn.dot(u, nrm)
    assert abs(s) < 0.5*eps, f"face is not planar s = {s}"
  return nrm, bar
  
def write_triangles_OBJ(wro, vtot, pock):
  # Writes the triangles of triangulations of the faces of the object to
  # OBJ file {wro}, assuming that the vertex indices are in the dict
  # {vtot}.
  
  assert not pock, "triangulaton with pockmarks not supported yet"
  
  # Returns a list {D} of all edges of all triangles, each 
  # as a pair {(kv_org,kv_dst)} with {kv_org < kv_dst}.
  # Element {D[0]} is not used.
  
  Nt_exp =  68 + 100 +  80   # Expected triangle count.
  Nd_exp = 3*Nt_exp // 2     # Expected triangulation edge count. 
  
  # Variables that are global to the nested functions:
  Ds = set()                # Triangulation edges.
  Tnrm = [ None ]           # Triangle normals.
  Tbar = [ None ]           # Triangle centers.
  debug_face = False 
  
  def Sved(p_org, p_dst):
    # Saves the edge from {p_org} to {p_dst} in {Ds},
    # if traversed in increasing vertex index sense.

    nonlocal Ds, Tnrm, Tbar, debug_face
    kv_org = p_org[3]
    kv_dst = p_dst[3]
    arc = (kv_org, kv_dst)
    if arc in Ds:
      sys.stderr.write(f"    !! duplicate arc {p_org[4]}--{p_dst[4]}\n")
      assert False
    Ds.add(arc)
    return None

  def Ptit(tit):
    # Title for a face or set of faces.

    nonlocal Ds, Tnrm, Tbar, debug_face
    wro.write("\n");
    wro.write("# %s\n" % tit);
    sys.stderr.write("writing %s \n" % (tit))
    return None
  
  def Prtt(key0, iv0, key1, iv1, key2, iv2, order):
    # Outputs triangle with vertices {vtot[key][iv]}
    # where {(key,iv)} is {(key0,iv0)}, {(key1,iv1)}, {(key2,iv2)}.
    # The vertices are listed in that order or the reverse
    # depending on whether {order} is {+1} or {-1}.

    nonlocal Ds, Tnrm, Tbar, debug_face
    # Save vertices in {Fp[0..2]} augmented with labels:
    Fp = []
    Fp.append(vtot[key0][iv0] + ( f"{key0}:{iv0:02d}", ))
    Fp.append(vtot[key1][iv1] + ( f"{key1}:{iv1:02d}", ))
    Fp.append(vtot[key2][iv2] + ( f"{key2}:{iv2:02d}", ))

    wro.write("f")
    for rv in range(3):
      jv = rv if order == +1 else 2-rv
      p = Fp[jv] 
      assert len(p) == 5
      kv = p[3]
      lab = p[4]
      if debug_face:
        sys.stderr.write(" v%04d = %s = ( %9.*f %9.*f %9.*f )\n" % (kv, lab, prec, p[0], prec, p[1], prec, p[2]));
      wro.write(" %d" % kv);
      Sved(Fp[jv], Fp[(jv + 3 + order) % 3])
    wro.write("\n");
    if debug_face: sys.stderr.write("\n")

    nrm, bar = face_normal_and_barycenter(Fp)
    Tnrm.append(nrm)
    Tbar.append(bar)
    return None
     
  # Here are the faces:
      
  Ptit(f"faces with normal ±X")
    
  for sx in +1, -1:
    tx = 'mop'[sx+1]

    Ptit(f" faces in halfspace sx = {sx:+d}")
    
    # Vert face of cube:
    Prtt(tx+'pm', 4,  tx+'pp', 4, tx+'pp', 7, sx)
    Prtt(tx+'mm', 4,  tx+'mm', 7, tx+'mp', 4, sx)
    Prtt(tx+'mm', 4,  tx+'pm', 4, tx+'pm', 7, sx)
    Prtt(tx+'mp', 7,  tx+'pp', 4, tx+'mp', 4, sx)
    Prtt(tx+'pp', 7,  tx+'pp', 4, tx+'mp', 7, sx)
    Prtt(tx+'pm', 4,  tx+'pp', 7, tx+'pm', 7, sx)
    Prtt(tx+'mm', 4,  tx+'pm', 7, tx+'mm', 7, sx)
    Prtt(tx+'mm', 7,  tx+'mp', 7, tx+'mp', 4, sx)
   
    # Centrad vert wall of hole on cube:
    Prtt(tx+'pp', 5,  tx+'mp', 5, tx+'mm', 5, sx)
    Prtt(tx+'pp', 5,  tx+'mm', 5, tx+'pm', 5, sx)
    
    # Distad vert wall of hole on cube:
    Prtt(tx+'mp', 6,  tx+'pp', 6, tx+'pm', 6, sx)
    Prtt(tx+'mp', 6,  tx+'pm', 6, tx+'mm', 6, sx)
    
    # Centrad vert wall of hole on handle:
    Prtt(tx+'pp', 9,  tx+'mp', 9, tx+'pm', 9, sx)
    Prtt(tx+'mm', 9,  tx+'pm', 9, tx+'mp', 9, sx)
    
    # Distad vert wall of hole on handle:
    Prtt(tx+'pp',10,  tx+'pm',10, tx+'mp',10, sx)
    Prtt(tx+'mp',10,  tx+'pm',10, tx+'mm',10, sx)
    
    # Extreme vert wall of handle:
    Prtt(tx+'pp', 8,  tx+'mp', 8, tx+'mm', 8, sx)
    Prtt(tx+'mm', 8,  tx+'pm', 8, tx+'pp', 8, sx)

    for sz in +1, -1:
      tz = 'mop'[sz+1]
    
      Ptit(f"    faces in quadrant sx = {sx:+d} sz = {sz:+d}")
      
      # Centrad vert wall of slot on face of cube:
      Prtt(tx+'m'+tz, 1,  tx+'p'+tz, 1,  tx+'m'+tz, 0, sx*sz)
      Prtt(tx+'p'+tz, 1,  tx+'p'+tz, 0,  tx+'m'+tz, 0, sx*sz)
      
      # Distad vert wall of slot on face of cube:
      debug_face = sx == -1 and sz == -1
      Prtt(tx+'p'+tz, 2,  tx+'m'+tz, 2,  tx+'m'+tz, 3, sx*sz)
      debug_face = False
      Prtt(tx+'p'+tz, 3,  tx+'p'+tz, 2,  tx+'m'+tz, 3, sx*sz)
      
      # Distad vert wall of bar in hole:
      Prtt(tx+'p'+tz,12,  tx+'m'+tz,12,  tx+'p'+tz,14, sx*sz)
      Prtt(tx+'p'+tz,14,  tx+'m'+tz,12,  tx+'m'+tz,14, sx*sz)
       
      # Centrad vert wall of bar in hole:
      Prtt(tx+'m'+tz,11,  tx+'p'+tz,11,  tx+'p'+tz,13, sx*sz)
      Prtt(tx+'m'+tz,13,  tx+'m'+tz,11,  tx+'p'+tz,13, sx*sz)
   
    
  Ptit(f"faces with normal ±Y")

  for sy in +1, -1:
    ty = 'mop'[sy+1]
    
    Ptit(f"  faces in halfspace sy = {sy:+d}")
    
    # Vert face of main cube:
    Ptit(f"    big complicated face")

    Prtt('m'+ty+'m', 4,  'm'+ty+'m', 2,  'm'+ty+'m', 3, sy)
    Prtt('m'+ty+'m', 4,  'm'+ty+'m', 6,  'm'+ty+'m', 2, sy)
    Prtt('m'+ty+'m', 4,  'm'+ty+'p', 6,  'm'+ty+'m', 6, sy)
    Prtt('m'+ty+'m', 4,  'm'+ty+'p', 4,  'm'+ty+'p', 6, sy)
    Prtt('m'+ty+'p', 6,  'm'+ty+'p', 4,  'm'+ty+'p', 2, sy)
    Prtt('m'+ty+'p', 2,  'm'+ty+'p', 4,  'm'+ty+'p', 3, sy)
                                      
    Prtt('m'+ty+'m', 2,  'm'+ty+'m', 5,  'm'+ty+'m', 1, sy)
    Prtt('m'+ty+'m', 2,  'm'+ty+'m', 6,  'm'+ty+'m', 5, sy)
    Prtt('m'+ty+'p', 6,  'm'+ty+'p', 2,  'm'+ty+'p', 5, sy)
    Prtt('m'+ty+'p', 5,  'm'+ty+'p', 2,  'm'+ty+'p', 1, sy)
                                      
    Prtt('m'+ty+'m', 0,  'm'+ty+'m', 1,  'o'+ty+'m', 0, sy)
    Prtt('o'+ty+'m', 0,  'm'+ty+'m', 1,  'o'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 1,  'm'+ty+'m', 1,  'm'+ty+'m', 5, sy)
    Prtt('o'+ty+'m', 1,  'm'+ty+'m', 5,  'm'+ty+'p', 5, sy)
    Prtt('o'+ty+'p', 0,  'm'+ty+'p', 5,  'm'+ty+'p', 1, sy)
    Prtt('o'+ty+'p', 0,  'm'+ty+'p', 1,  'm'+ty+'p', 0, sy)
                                      
    Prtt('m'+ty+'p', 5,  'o'+ty+'p', 0,  'p'+ty+'p', 5, sy)
    Prtt('p'+ty+'p', 5,  'o'+ty+'m', 1,  'm'+ty+'p', 5, sy)
                                      
    Prtt('p'+ty+'m', 0,  'o'+ty+'m', 0,  'p'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 0,  'o'+ty+'m', 1,  'p'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 1,  'p'+ty+'m', 5,  'p'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 1,  'p'+ty+'p', 5,  'p'+ty+'m', 5, sy)
    Prtt('o'+ty+'p', 0,  'p'+ty+'p', 1,  'p'+ty+'p', 5, sy)
    Prtt('o'+ty+'p', 0,  'p'+ty+'p', 0,  'p'+ty+'p', 1, sy)
                                      
    Prtt('p'+ty+'m', 1,  'p'+ty+'m', 5,  'p'+ty+'m', 2, sy)
    Prtt('p'+ty+'m', 2,  'p'+ty+'m', 5,  'p'+ty+'m', 6, sy)
    Prtt('p'+ty+'p', 5,  'p'+ty+'p', 2,  'p'+ty+'p', 6, sy)
    Prtt('p'+ty+'p', 2,  'p'+ty+'p', 5,  'p'+ty+'p', 1, sy)
                                      
    Prtt('p'+ty+'m', 4,  'p'+ty+'m', 3,  'p'+ty+'m', 2, sy)
    Prtt('p'+ty+'m', 4,  'p'+ty+'m', 2,  'p'+ty+'m', 6, sy)
    Prtt('p'+ty+'m', 4,  'p'+ty+'m', 6,  'p'+ty+'p', 6, sy)
    Prtt('p'+ty+'m', 4,  'p'+ty+'p', 6,  'p'+ty+'p', 4, sy)
    Prtt('p'+ty+'p', 4,  'p'+ty+'p', 6,  'p'+ty+'p', 2, sy)
    Prtt('p'+ty+'p', 4,  'p'+ty+'p', 2,  'p'+ty+'p', 3, sy)
    
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sx = {sx:+d} sy = {sy:+d}")
     
      # Vert outer face of handle:
      Prtt(tx+ty+'p', 7,  tx+ty+'p', 8,  tx+ty+'m', 8, sx*sy)
      Prtt(tx+ty+'p', 7,  tx+ty+'m', 8,  tx+ty+'m', 7, sx*sy)
                                      
      # Vert inner face of handle:
      Prtt(tx+ty+'p',10,  tx+ty+'p', 9,  tx+ty+'m', 9, sx*sy)
      Prtt(tx+ty+'p',10,  tx+ty+'m', 9,  tx+ty+'m',10, sx*sy)

      for sz in +1, -1:
        tz = 'mop'[sz+1]

        Ptit(f"      faces in octant sx = {sx:+d} sy = {sy:+d} sz = {sz:+d}")

        # Tip of bar in hole:
        Prtt(tx+ty+tz,13,  tx+ty+tz,11,  tx+ty+tz,14, sx*sy*sz)
        Prtt(tx+ty+tz,14,  tx+ty+tz,11,  tx+ty+tz,12, sx*sy*sz)
      
    
  Ptit(f"faces with normal ±Z")

  for sz in +1, -1:
    tz = 'mop'[sz+1]
    
    Ptit(f"  faces in halfspace sz = {sz:+d}")
    
    # Centrad horz face of main cube:
    Prtt('op'+tz, 0,  'om'+tz, 0,  'pp'+tz, 0, sz)
    Prtt('om'+tz, 0,  'pm'+tz, 0,  'pp'+tz, 0, sz)
    Prtt('om'+tz, 0,  'op'+tz, 0,  'mm'+tz, 0, sz)
    Prtt('op'+tz, 0,  'mp'+tz, 0,  'mm'+tz, 0, sz)
    
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sx = {sx:+d} sz = {sz:+d}")

      # Distad horz face of main cube:
      Prtt(tx+'m'+tz, 3,  tx+'m'+tz, 4,  tx+'p'+tz, 4, sx*sz)
      Prtt(tx+'m'+tz, 3,  tx+'p'+tz, 4,  tx+'p'+tz, 3, sx*sz)
      
      # Horz floor face of slot on face of cube:
      Prtt(tx+'m'+tz, 1,  tx+'m'+tz, 2,  tx+'p'+tz, 2, sx*sz)
      Prtt(tx+'m'+tz, 1,  tx+'p'+tz, 2,  tx+'p'+tz, 1, sx*sz)
      
      # Distad floor of bar in hole:
      Prtt(tx+'p'+tz,12,  tx+'p'+tz,11,  tx+'m'+tz,12, sx*sz)
      Prtt(tx+'p'+tz,11,  tx+'m'+tz,11,  tx+'m'+tz,12, sx*sz)
       
      # Centrad floor of bar in hole:
      Prtt(tx+'p'+tz,13,  tx+'p'+tz,14,  tx+'m'+tz,13, sx*sz)
      Prtt(tx+'p'+tz,14,  tx+'m'+tz,14,  tx+'m'+tz,13, sx*sz)

      # Horz face of hole on main cube:
      Prtt(tx+'p'+tz, 5,  tx+'p'+tz, 6,  tx+'m'+tz, 5, sx*sz)
      Prtt(tx+'m'+tz, 5,  tx+'p'+tz, 6,  tx+'m'+tz, 6, sx*sz)
      
      # Horz face on handle:
      Prtt(tx+'p'+tz, 9,  tx+'p'+tz, 8,  tx+'p'+tz, 7, sx*sz)
      Prtt(tx+'p'+tz, 9,  tx+'p'+tz,10,  tx+'p'+tz, 8, sx*sz)
      Prtt(tx+'m'+tz, 8,  tx+'p'+tz, 8,  tx+'p'+tz,10, sx*sz)
      Prtt(tx+'m'+tz, 8,  tx+'p'+tz,10,  tx+'m'+tz,10, sx*sz)
      Prtt(tx+'m'+tz, 7,  tx+'m'+tz, 8,  tx+'m'+tz,10, sx*sz)
      Prtt(tx+'m'+tz, 7,  tx+'m'+tz,10,  tx+'m'+tz, 9, sx*sz)
      Prtt(tx+'m'+tz, 7,  tx+'m'+tz, 9,  tx+'p'+tz, 7, sx*sz)
      Prtt(tx+'m'+tz, 9,  tx+'p'+tz, 9,  tx+'p'+tz, 7, sx*sz)
    
  Nt_cmp = len(Tnrm) - 1
  sys.stderr.write("wrote %d triangles (expected %d)\n" % (Nt_cmp, Nt_exp))
  assert Nt_cmp == Nt_exp
    
  # Convert to list, leaving only one arc per edge:
  D = [ None ] + [ arc for arc in Ds if arc[0] < arc[1] ]
  
  Nd_cmp = len(D) - 1
  sys.stderr.write("wrote %d edges (expected %d)\n" % (Nd_cmp, Nd_exp))
  #assert Nd_cmp == Nd_exp

  return Tnrm, Tbar, D
  
def face_normal_and_barycenter(Fp):
  # Compute normal {nrm} and barycenter {bar}
  # given a list of vertices {Fp[0..deg-1]}
  deg = len(Fp)
  assert deg >= 3
  po = Fp[0][:3]
  pa = Fp[1][:3]
  pb = Fp[2][:3]
  u = rn.sub(pa,po)
  v = rn.sub(pb,po)
  nrm, sz = rn.dir(rn.cross3d(u, v))

  bar = (0, 0, 0)
  for jv in range(deg):
    pj = Fp[jv][:3]
    bar = rn.mix(1.0, bar, 1.0/deg, pj)
    
  # Check planarity:
  for jv in range(deg):
    pj = Fp[jv][:3]
    u = rn.sub(pj, bar)
    s = rn.dot(u, nrm)
    assert abs(s) < 0.5*eps, f"face is not planar s = {s}"
  return nrm, bar

def rotation_matrix(xrot,yrot):
  # Returns a 3x3 orthonormal matrix that
  # rotates by {xrot} degrees towatds the {X}-axis
  # and {yrot} degrees toward the {Y}-axis.
  
  if xrot == 0 and yrot == 0:
    rmat = rmxn.ident_matrix(3,3)
  else:
    t, et = rn.dir((tan(xrot*pi/180), tan(yrot*pi/180), 1))
    assert et != 0
    s, es = rn.dir(rn.cross3d((0,0,1), t))
    assert es != 0
    r, er = rn.dir(rn.cross3d(t, s))
    assert abs(er - 1.0) <= 0.001
    rmat = (r, s, t)
    
  return rmat

def reven(C):
  # Rounds {C} to even multiple of {eps}.
  return eps*2*floor(C/eps/2 + 0.5)

def radians(deg):
  # Converts degrees to radians.
  return pi*deg/180

def checkpow(n, msg):
  # Checks that {n} is a power of 2:
  pp = n;
  while pp % 2 == 0: pp = pp//2
  assert pp == 1, msg

def main():
  pert = int(sys.argv[1])          # Random perturbation amount, in {eps} units.
  xrot = int(sys.argv[2])          # Rotation in the direction of the {X}-axis (degrees).
  yrot = int(sys.argv[3])          # Rotation in the direction of the {Y}-axis (degrees).
  pock = (int(sys.argv[4]) != 0)   # Add pockmarks to big complex face.
  
  sys.stderr.write("slicing_hel_example.py: ")
  sys.stderr.write(f" pert = {pert}")
  sys.stderr.write(f" xrot = {xrot:.4f}")
  sys.stderr.write(f" yrot = {yrot:.4f}")
  sys.stderr.write(f" pock = {pock}")
  sys.stderr.write("\n")
  
  assert 0 <= pert and pert < 50, "invalid pert"
  assert abs(xrot) < 90, "invalid xrot"
  assert abs(yrot) < 90, "invalid yrot"

  S = 100       # Side of cube
  
  rmat = rotation_matrix(xrot,yrot)

  vref = ref_vertices(S, pock)
  vtot, Nv = all_vertices(vref)
  perturb_and_round_vertices(vtot, rmat,pert)

  file_prefix = f"out/hel_pa{pert:03d}_xr{xrot:03d}_yr{yrot:03d}_pk{str(pock)[0]}"

  wro = open(file_prefix + "_triF.obj", 'w')
  write_vertices_OBJ(wro, vtot)
  Fnrm, Fbar, E, G = write_faces_OBJ(wro, vtot, pock)
  Ne = len(E) - 1
  Ng = len(G) - 1
  Nf = len(Fnrm) - 1
  wro.close()

  if not pock:
    wrt = open(file_prefix + "_triT.obj", 'w')
    write_vertices_OBJ(wrt, vtot)
    Tbar, Tnrm, D = write_triangles_OBJ(wrt, vtot, pock)
    Nd = len(D) - 1
    Nt = len(Tnrm) - 1
    wrt.close()

    # Exclude real and ghost edges from the {D} list:
    EGs = set((E[1:] + G[1:]))
    Ds = set(D[1:])
    Ds = Ds - EGs
    D = [None] + list(Ds)
    Nd = len(D) - 1 

    write_povray = False

    if write_povray:
      wrp = open(file_prefix + ".inc", 'w')
      povray.write_parms(wrp, S, Nv,Ne,Ng,Nf,Nt,Nd)
      povray.write_vertices(wrp, vobj, Vlab)
      povray.write_edges(wrp, E, G, D, vobj, Vlab)
      povray.write_faces(wrp, Fnrm, Fbar)
      wrp.close()
  
  return 0

main()