#! /usr/bin/python3
# Last edited on 2024-05-30 09:03:21 by stolfi

from math import sin, cos, tan, sqrt, pi, inf, floor
import rn
import sys
import re
 
prec = 4          # Decimal digits    
eps = 0.1**prec   # Coordinate quantum.

# The object is a nut-like solid that is basically the difference
# between two concentric prisms with {Nr} sides, with chamfered edges.
# 
# The set of vertices consists of {2*Nr} trozos indexed by indices {kr}
# in {0..Nr-1} and index {km} in {0,1}. Trozo {[kr][km]} differs from
# trozo {0,0} by a rotation of {2*pi*kr/Nr} radians about {Z} and
# mirroring {km} times about the plane {Z = 0}.
#
# Each trozo has an outer corner vertex {va}, {Nc+1} lateral chamfer
# vertices {vw[0..Nc]}, {Nc+1} top chamfer vertices {vt[0..Nc]}, and a
# top hole vertex {vh}.

def ref_trozo_vertices(Nr, Nc, Ro, Ri, Rc, Ha):
  # Computes the vertices of trozo {[0][0]}. Parameters:
  
  # {Nr} Num sides of nut.
  # {Nc} Num facets in chamfer.
  # {Ro} Circumradius of nut.
  # {Ri} Circumradius of nut hole.
  # {Rc} Radius of chamfer vertices on top face.
  # {Ha} Total actual height of nut.
  
  # These vertices are returned as a {dict} with keys {'va','vw','vt','vh'}.
  
  # Key derived dimensions (before from chamfering):
  hang = pi/Nr           # Half the angular span of each trozo.
  Xoface = Ro*cos(hang)  # {X} coordinate of outer wall.
  Yoface = Ro*sin(hang)  # Max abs {Y} of outer wall.

  Ztop = Ha/2             # {Z} coordinate of nut's top face.

  # The /outer wall/ is the flat outer face of trozo. The /ideal outer
  # wall/ would be a rectangle on the plane {X = Xoface} extending in {Y}
  # over {[-Yoface _ +Yoface]} and in {Z} over {[0 _ Ztop]}. The /outer
  # ridge/ is the upper side of this rectangle, where it meets the top
  # plane of the nut. The /outer quines/ are the vertical sides of the
  # rectangle, where the ideal outer wall of this trozo meets the ideal
  # outer walls of adjacent trozos.

  Xiface = Ri*cos(hang)  # {X} coordinate of inner wall.
  Yiface = Ri*sin(hang)  # Max abs {Y} of inner wall.

  # The /inner wall/ is the vertical face of the trozo that is one of 
  # the {Nr} faces of the hole. Since there is no chamfer there,
  # it would be a rectangle on the plane {X = Xiface} extending in {Y}
  # over {[-Yiface _ +Yiface]} and in {Z} over {[0 _ Ztop]}. The 
  # /inner ridge/ is the side of the rectangle at {Z = Ztop},
  # and the /inner quines/ are the sides parallel to the {Z} axis. 

  sys.stderr.write("Nr = %d Nc = %d\n" % (Nr, Nc))
  sys.stderr.write("Ro = %.*f Ri = %.*f Ha = %.*f\n" % (prec, Ro, prec, Ri, prec, Ha))
  
  # The /actual outer wall/ of the trozo is what remains of the 
  # ideal outer wall after chamfering is applied.  The chamfering eats away
  # the whole outer ridge and the top part of the outer quines.
  # The chamfering is a band of {Nc} quadrangular facets that bridges 
  # two pentagonal /corner facets/ shared with adjacent trozos.

  # Compute chamfer vertices {vw[kc],vt[kc]}:
  dach = 2*hang/(Nc+1)          # Angle between chamfer vertices on top face.
  vw = [ None ]*(Nc+1)
  vt = [ None ]*(Nc+1)
  for kc in range(Nc+1):
    ach = (kc - Nc/2)*dach       # Argument of chamfer vertex on top face.
    Yuv = Rc*sin(ach)            # {Y} coordinate of {vw[kc],vt[kc]}.
    Xu = Rc*cos(ach)             # {X} coord of {vt[kc]}.
    dXZ = Xoface - Xu            # {X} distance from top chamfer vertex and wall.  
    Zv = Ztop - dXZ              # {Z} coord of {vw[kc]}
    
    vt[kc] = ( Xu, Yuv, Ztop )
    vw[kc] = ( Xoface, Yuv, Zv )
    
  # Compute the coeffs of {Z = A*X + B*Y + C} for plane of pentagonal facet:
  Xchmax = vt[0][0]         # {X} coordinate of first chamfer vertex on top face.
  Ychmax = vt[0][1]         # {Y} coordinate of first chamfer vertex on top face.
  dXZmax = Xoface - Xchmax  # Max chamfer vertex {X},{Z} distance from outer ridge.

  A = -1
  B = tan(hang)
  C = Ztop - A*Xchmax - B*Ychmax;
  # Check equation:
  Zcheck = A*Xoface + B*Ychmax + C
  assert abs(Zcheck - (Ztop - dXZmax)) < 0.0001, "bad math"

  # Compute corner vertex from equation:
  Xa = Xoface
  Ya = -Yoface
  Za = A*Xa + B*Ya + C
  va = ( Xa, Ya, Za )
  
  # Compute the top inner hole corner vertex:
  Yh = -Yiface
  vh = ( Xiface, Yh, Ztop )
  
  Tref = { 'va': va, 'vh': vh, 'vw': vw, 'vt': vt }
  return Tref
  
def all_vertices(Tref, Nr):
  # The structure {Tref} must be a {dict} that gives the vertices of the
  # reference trozo of the object, as produced by {ref_trozo_vertices}.
  # Returns an array {Stot[0..Nr-1][0..1]} where each element is a copy
  # of {Tref} properly rotated and mirrored in {Z}.

  ang = 2*pi/Nr
  sa = sin(ang); ca = cos(ang)
  
  def rotate_point(p):
    return ( ca*p[0] - sa*p[1], sa*p[0] + ca*p[1], p[2] )

  def mirror_point(p):
    return ( p[0], p[1], -p[2] )

  Stot = alloc_array(Nr, 2)
  Trot = Tref;
  for kr in range(Nr):
    assert isinstance(Trot, dict)
    Trot = map_trozo(Trot, rotate_point)
    Stot[kr][0] = Trot
    Stot[kr][1] = map_trozo(Trot, mirror_point)
  return Stot
  
def alloc_array(m, n):
  # Returns a list of {m} lists of {n} elements, all {None}.
  A = [ None ]*m
  for i in range(m):
    A[i] = [ None]*n
  return A
  
def make_null_trozo(Nc):
  # Creates an array like the total vertex table but with all
  # entries set to {None}.
  Vnull = {}
  Vnull['va'] = None
  Vnull['vh'] = None
  Vnull['vw'] = [ None ]*(Nc+1)
  Vnull['vt'] = [ None ]*(Nc+1)
  return Vnull

def map_trozo(V, func):
  Nc = len(V['vw']) - 1;
  Vmap = make_null_trozo(Nc)
  Vmap['va'] = func(V['va'])
  Vmap['vh'] = func(V['vh'])
  for kc in range(Nc+1):
    Vmap['vw'][kc] = func(V['vw'][kc])
    Vmap['vt'][kc] = func(V['vt'][kc])
  return Vmap
  
def write_vertices_OBJ(wro, Stot, Nr):
  # Expect {Stot} to be a an array of {Nr × 2} trozos comprising the
  # vertices of a nut object. Writes the vertices to {wro} in OBJ
  # format.
  # 
  # Rounds each coordinate to even multiple of {eps}. Returns an array
  # {Sind} with the same shape as {Stot} but where each point is
  # replaced by the integer index of the vertex, as expected by the OBJ
  # format (from 1).
  #
  # Also returns a list {Vpos} where element {Vpos[i]} is the
  # coordinates of the vertex with OBJ index {i}, for {i} in {1..Nv}
  # where {Nv} is the total number of vertices. Note that {Vpos[0]} is
  # undefined.
  # 
  # Also returns a list {Vlab} where element {Vlab[i]} is the
  # position in {Stot} of vertex {Vpos[i]}. For example "3.0.va" to mean
  # vertex {'va'} of {Stot[3][0]['va']}, and 
  # "3.0.vt.5" to mean {Stot[3][0]['vt'][5]}.

  Nc = len(Stot[0][0]['vw']) - 1  # Number of chamfer facets per trozo.
  Nv = 2*Nr*(2*(Nc+1) + 2)        # Expected number of vertices in object.
  
  Vpos = [ None ] # Vertex positions. Element 0 is not used.
  Vlab = [ None ] # Vertex labels. # Element 0 is not used.
    
  def Prtv(p,lab):
    # Writes {p} to {wro}, with coords rounded to even multiples of {eps}. 
    # Also appends it to {Vpos,Vlab} (without rounding).
    nonlocal Vpos, Vlab
    wro.write("v");
    for i in range(3): wro.write(" %.*f" % (prec, reven(p[i])))
    wro.write("\n");
    Vpos.append(p)
    Vlab.append(lab)
    
  Sind = alloc_array(Nr, 2)

  kv = 0;
  for kr in range(Nr):
    for km in range(2):
      Sind[kr][km] = make_null_trozo(Nc)
      kv += 1;
      Prtv(Stot[kr][km]['va'], f"{kr}.{km}.va")
      Sind[kr][km]['va'] = kv
      kv += 1;
      Prtv(Stot[kr][km]['vh'], f"{kr}.{km}.vh")
      Sind[kr][km]['vh'] = kv
      for kc in range(Nc+1):
        kv += 1;
        Prtv(Stot[kr][km]['vw'][kc], f"{kr}.{km}.vw.{kc}")
        Sind[kr][km]['vw'][kc] = kv
        kv += 1;
        Prtv(Stot[kr][km]['vt'][kc], f"{kr}.{km}.vt.{kc}")
        Sind[kr][km]['vt'][kc] = kv

  assert kv == Nv
  assert len(Vpos) == Nv + 1
  assert len(Vlab) == Nv + 1

  wro.write("\n")

  return Sind, Vpos, Vlab
  
def write_faces_OBJ(wro, Sind, Vpos, Nr):
  # Writes the faces of the nut object {to OBJ file {wro}, assuming that the vertex indices are
  # in {Sind}. 
  
  # Also returns tables {Fnrm,Fpos} where {Fnrm[kf]} is the outward normal
  # of face {kf} and {Fpos[kf]} is a point on it, for {kf} in {1..Nf}.
  #
  # Also returns tables {Eorg,Edst} where {Eorg[ke]} and {Edst[ke]} are the
  # OBJ indices of the vertices which are the endpoints of edge {ke},
  # for {ke} in {1..Ne}.
  
  Nc = len(Sind[0][0]['vw']) - 1; # Number of chamfer facets per trozo.
  Nv = len(Vpos) - 1              # Total number of vertices.
  Nf = Nr*(8 + 2*Nc)              # Expected total number of face planes.
  Ne = Nr*(2*(3*Nc+7) + 2)        # Expected total number of edges.
  
  Fnrm = [ None ] # Outward normals of face planes, {[1..Nf]}.
  Fpos = [ None ] # A point on each face plane, {[1..Nf]}.
  Eorg = [ None ] # Indices of origin vertices of the edges, {[1..Ne]}.
  Edst = [ None ] # Indices of destination vertices of the edges, {[1..Ne]}.

  # Variables that are global to the nested functions:
  debug_face = False
  Fv = []    # Indices of vertices of current face.
  Fn = None  # Normal of current face.

  def Ptit(tit, kr, debug = False):
    # Title for a face or set of faces.
    nonlocal debug_face, Fv, Fn
    wro.write("\n");
    wro.write("# %s kr = %d\n" % (tit, kr));
    sys.stderr.write("# %s kr = %d\n" % (tit, kr));
    debug_face = debug
  
  def Bof():
    #  Start of a new face.
    nonlocal debug_face, Fv, Fn
    wro.write("f")
    assert len(Fv) == 0
    assert Fn == None
  
  def Prtv(iv):
    # Adds vertex with OBJ index {iv} to the current face.
    nonlocal debug_face, Fv, Fn
    wro.write(" %d" % iv);
    p = Vpos[iv]
    if debug_face:
      sys.stderr.write(" v%4d = ( %9.*f %9.*f %9.*f )\n" % (iv, prec, p[0], prec, p[1], prec, p[2]));
    # Saves vertex in list {Fv} of face vertices:
    Fv.append(iv);
    if len(Fv) == 3:
      # Compute normal {Fn}, appends to {Fnrm,Fpos}:
      po = Vpos[Fv[0]]
      pa = Vpos[Fv[1]]
      pb = Vpos[Fv[2]]
      u = rn.sub(pa,po)
      v = rn.sub(pb,po)
      Fn, sz = rn.dir(rn.cross3d(u, v))
      Fnrm.append(Fn)
      Fpos.append(po)
    elif len(Fv) > 3:
      # Check planarity:
      u = rn.sub(p, Vpos[Fv[0]])
      s = rn.dot(u, Fn)
      assert abs(s) < 0.5*eps, f"face is not planar s = {s}"
    if len(Fv) >= 2 and Fv[-2] < Fv[-1]:
      # Saves edge endpoints in {Eorg,Edst}:
      Eorg.append(Fv[-2])
      Edst.append(Fv[-1])
    
  def Eof():
    # End of face.
    nonlocal debug_face, Fv, Fn
    wro.write("\n");
    assert len(Fv) >= 3
    if Fv[-1] < Fv[0]:
      Eorg.append(Fv[-1])
      Edst.append(Fv[0])
    Fv = []; Fn = None
  
  for kr in range(Nr):
    krp = (kr + 1) % Nr # Index of next trozo pairs.
    krm = (kr - 1) % Nr # Index of previous trozo pair.
    assert krm >= 0 and krm < Nr
    
    Ptit("Outer wall face:", kr)
    Bof()
    Prtv(Sind[kr][0]['va'])
    Prtv(Sind[kr][1]['va'])
    for kc in range(Nc+1):
      Prtv(Sind[kr][1]['vw'][kc])
    Prtv(Sind[krp][1]['va'])
    Prtv(Sind[krp][0]['va'])
    for kc in range(Nc+1):
      Prtv(Sind[kr][0]['vw'][Nc-kc])
    Eof()
    
    Ptit("Top corner chamfer facet:", kr)
    Bof()
    Prtv(Sind[krm][0]['vt'][Nc])
    Prtv(Sind[krm][0]['vw'][Nc])
    Prtv(Sind[kr][0]['va'])    
    Prtv(Sind[kr][0]['vw'][0])
    Prtv(Sind[kr][0]['vt'][0])
    Eof()
    
    Ptit("Bottom corner chamfer facet:", kr)
    Bof()
    Prtv(Sind[krm][1]['vw'][Nc])
    Prtv(Sind[krm][1]['vt'][Nc])
    Prtv(Sind[kr][1]['vt'][0])
    Prtv(Sind[kr][1]['vw'][0])
    Prtv(Sind[kr][1]['va'])
    Eof()
    
    Ptit("Top chamfer facets:", kr)
    for kc in range(Nc):
      Bof()
      Prtv(Sind[kr][0]['vt'][kc])
      Prtv(Sind[kr][0]['vw'][kc])
      Prtv(Sind[kr][0]['vw'][kc+1])
      Prtv(Sind[kr][0]['vt'][kc+1])
      Eof()
    
    Ptit("Bottom chamfer facets:", kr)
    for kc in range(Nc):
      Bof()
      Prtv(Sind[kr][1]['vt'][kc+1])
      Prtv(Sind[kr][1]['vw'][kc+1])
      Prtv(Sind[kr][1]['vw'][kc])
      Prtv(Sind[kr][1]['vt'][kc])
      Eof()
      
    Ptit("Top face between chamfer and hole:", kr)
    Bof()
    Prtv(Sind[kr][0]['vh'])
    for kc in range(Nc+1):
      Prtv(Sind[kr][0]['vt'][kc])
    Prtv(Sind[krp][0]['vh'])
    Eof()
      
    Ptit("Bottom face between chamfer and hole:", kr)
    Bof()
    Prtv(Sind[krp][1]['vh'])
    for kc in range(Nc+1):
      Prtv(Sind[kr][1]['vt'][Nc-kc])
    Prtv(Sind[kr][1]['vh'])
    Eof()
      
    Ptit("Top face between corner chamfer face and hole:", kr)
    Bof()
    Prtv(Sind[krm][0]['vt'][Nc])
    Prtv(Sind[kr][0]['vt'][0])
    Prtv(Sind[kr][0]['vh'])
    Eof()
      
    Ptit("Bottom face between corner chamfer face and hole:", kr)
    Bof()
    Prtv(Sind[kr][1]['vh'])
    Prtv(Sind[kr][1]['vt'][0])
    Prtv(Sind[krm][1]['vt'][Nc])
    Eof()
    
  # Hole wall faces must be last, for POV-Ray's sake:
  for kr in range(Nr):
    krp = (kr + 1) % Nr # Index of next trozo pairs.
    krm = (kr - 1) % Nr # Index of previous trozo pair.
    assert krm >= 0 and krm < Nr
    
    Ptit("Inner wall face:", kr)
    Bof()
    Prtv(Sind[kr][0]['vh'])
    Prtv(Sind[kr][1]['vh'])
    Prtv(Sind[krm][1]['vh'])
    Prtv(Sind[krm][0]['vh'])
    Eof()
    
  sys.stderr.write("wrote %d faces (expected %d)\n" % (len(Fnrm)-1, Nf))
  assert len(Fnrm) == Nf + 1
  assert len(Fpos) == Nf + 1
  
  sys.stderr.write("found %d edges (expected %d)\n" % (len(Eorg)-1, Ne))
  assert len(Eorg) == Ne + 1
  assert len(Edst) == Ne + 1

  return Fnrm, Fpos, Eorg, Edst
  
def write_vertices_POV(wrp, Vpos, Vlab):
  # Writes to {wrp} a POV-Ray macro that defines 
  # the vertices {Vpos[1..Nv]} in POV-ray format.
  # The labels {Vlab[1..Nv]} are written as comments.

  Nv = len(Vpos) - 1
  
  # Write the POV-Ray preamble:
  wrp.write("#macro nut_vertices()\n")
  wrp.write("  #local Nv = nut_num_vertices;\n")
  wrp.write("  // Returns an array {V} of {Nv+1} elements.\n")
  wrp.write("  // Element {V[kv]} is the vertex with OBJ index {kv} in {1..Nv}.\n")
  wrp.write("  // Element {V[0]} is not used.\n")
  wrp.write("  #local V = array[Nv+1]\n")
  wrp.write("  #local V[   0] = < -1, -1, -1>; // Not used.\n")

  for kv in range(1,Nv+1):
    p = Vpos[kv]
    lab = Vlab[kv]
    wrp.write("  #local V[%4d] =" % kv)
    wrp.write(" < %.*f, %.*f, %.*f >;" % (prec, reven(p[0]), prec, reven(p[1]), prec, reven(p[2])))
    wrp.write(" // %s\n" % lab)

  # Write the POV-Ray postamble:
  wrp.write("  V\n")
  wrp.write("#end\n\n")
  
def write_edges_POV(wrp, Eorg, Edst, Nr, Vpos, Vlab):
  # Writes to {wrp} the indices of the endpoints of the edges of the nut.
  
  Nv = len(Vpos) - 1
  Ne = len(Eorg) - 1

  # Write the POV-Ray preambles:
  wrp.write("#macro nut_edges()\n")
  wrp.write("  #local Ne = nut_num_edges;\n")
  wrp.write("  // Returns an array {E} of {Ne+1} elements.\n")
  wrp.write("  // Element {E[ke]} is the triple {<ko[ke],kd[ke],ty[ke]>}, for {ke} in {1..Ne}.\n")
  wrp.write("  // Here {ko[ke],kd[ke]} are the indices of the vertices\n")
  wrp.write("  // which are the endpoints of edge {E[ke]},\n")
  wrp.write("  // and {ty[ke]} is an edge type code:\n")
  wrp.write("  // \n")
  wrp.write("  //   100-199 outer vertical edge.\n")
  wrp.write("  //   200-299 inner (hole) vertical edge.\n")
  wrp.write("  //   300-399 outer wall bottom edge.\n")
  wrp.write("  //   400-499 inner (hole) wall bottom edge.\n")
  wrp.write("  //   0       other.\n")
  wrp.write("  // \n")
  wrp.write("  // The rotation index of the edge is {ty[ke] % 100}. For quines,\n")
  wrp.write("  // the rotation index is ???.\n")
  wrp.write("  #local E = array[Ne+1]\n")
  wrp.write("  #local E[   0] = < -1, -1, -1 >; // Not used.\n")
  
  for ke in range(1, Ne+1):
    ko = Eorg[ke];
    kd = Edst[ke]; 
    
    # Determine edge type from labels of its endpoints:
    labo = Vlab[ko]; kro = int(re.sub(r"[.].*$", "", labo)) # Label and rotation index of {ko}.
    labd = Vlab[kd]; krd = int(re.sub(r"[.].*$", "", labd)) # Label and rotation index of {kd}.
    vw1_pat = "^[0-9]+[.]1[.]vw.*"    # Bottom outer wall chamfer vertex.
    vw1_o = re.match(vw1_pat, labo); vw1_d = re.match(vw1_pat, labd)

    vh1_pat = "^[0-9]+[.][1][.]vh.*"  # Vertex on bottom of hole wall.
    vh1_o = re.match(vh1_pat, labo); vh1_d = re.match(vh1_pat, labd)
    
    va_pat = "^[0-9]+[.][01][.]va.*"  # Vertex on vertical side of outer wall. 
    va_o = re.match(va_pat, labo); va_d = re.match(va_pat, labd)
    
    vh_pat = "^[0-9]+[.][01][.]vh.*"  # Vertex on vertical side of hole wall.
    vh_o = re.match(vh_pat, labo); vh_d = re.match(vh_pat, labd)

    if va_o and va_d and kro == krd:
      # Edge is an outer quine:
      ty = 100 + kro
    elif vh_o and vh_d and kro == krd:
      # Edge is an inner (hole) quine:
      ty = 200 + kro
    elif vw1_o and (vw1_d or va_d):
      # Edge is on lower border of outer wall:
      ty = 300 + kro
    elif vw1_d and (vw1_o or va_o):
      # Edge is on lower border of outer wall:
      ty = 300 + krd
    elif vh1_o and vh1_d and kro != krd:
      # Edge is on lower border of inner wall:
      if kro == (krd+1) % Nr:
        ty = 400 + krd
      elif krd == (kro+1) % Nr:
        ty = 400 + kro
    else:
      ty = 0

    wrp.write("  #local E[%4d] = < %d, %d, %d >;\n" % (ke, Eorg[ke], Edst[ke], ty))
    elen = rn.norm(rn.sub(Vpos[Edst[ke]], Vpos[Eorg[ke]]))
    assert elen >= 5*eps, "edge too short"
    
  # Write the POV-Ray postamble:
  wrp.write("  E\n")
  wrp.write("#end\n\n")
 
def write_faces_POV(wrp, Fnrm, Fpos, Nr):
  # Writes the face planes to {wrp} as POV-ray.
  
  Nf = len(Fnrm) - 1

  # Write the POV-Ray preambles:
  wrp.write("#macro nut_faces()\n")
  wrp.write("  #local Nf = nut_num_faces;\n")
  wrp.write("  #local Nr = nut_num_hole_sides;\n")
  wrp.write("  // Returns an array {F} of {Nf+1} elements.\n")
  wrp.write("  // Element {F[kf]} is an instance of the\n")
  wrp.write("  // 'plane' (halfspace) POV-ray primitive\n")
  wrp.write("  // for {kf} in {1..Nf}.\n")
  wrp.write("  // The last {Nr} planes are walls of the hole.\n")
  wrp.write("  #local F = array[Nf+1]\n")
  wrp.write("  #local F[   0] = sphere{ <0,0,0>, 1000} // Not used.\n")
  
  for kf in range(1, Nf+1):
    Fn = Fnrm[kf]
    Fp = Fpos[kf]
    wrp.write("  #local F[%4d] = plane{ " % kf)
    wrp.write(" < %.7f, %.7f, %.7f >, 0" % (Fn[0], Fn[1], Fn[2]))
    wrp.write(" translate < %.*f, %.*f, %.*f >" % (prec, Fp[0], prec, Fp[1], prec, Fp[2]))
    wrp.write(" }\n")
    
  # Write the POV-Ray postamble:
  wrp.write("  F\n")
  wrp.write("#end\n\n")

def main():
  Nr = 9                  # Num sides of nut.  Must be small.
  Nc = int(sys.argv[1])   # Num facets in each ridge chamfer.
  Ro = 100                # Circumradius of nut.
  Ri = 70                 # Circumradius of nut hole.
  Rc = 0.95*Ro*cos(pi/Nr) # Radius of chamfer vertices on top face.
  Hn = 40                 # Nominal height of nut.
  
  Ha = Hn - 4*eps         # Actual height of nut.

  Tref = ref_trozo_vertices(Nr, Nc, Ro, Ri, Rc, Ha)
  Stot = all_vertices(Tref, Nr)

  file_pref = f"out/nut_{Nc:08d}"

  wro = open(file_pref + ".obj", 'w')
  Sind, Vpos, Vlab = write_vertices_OBJ(wro, Stot, Nr)
  Fnrm, Fpos, Eorg, Edst = write_faces_OBJ(wro, Sind, Vpos, Nr)
  wro.close()
  
  Nv = len(Vpos) - 1
  Ne = len(Eorg) - 1
  Nf = len(Fnrm) - 1
  
  wrp = open(file_pref + ".inc", 'w')
  wrp.write("#declare nut_outer_radius = %.*f;\n\n" % (prec, Ro))
  wrp.write("#declare nut_hole_radius = %.*f;\n\n" % (prec, Ri))
  wrp.write("#declare nut_height = %.*f;\n\n" % (prec, Hn))
  wrp.write("#declare nut_num_faces = %d;\n\n" % Nf)
  wrp.write("#declare nut_num_edges = %d;\n\n" % Ne)
  wrp.write("#declare nut_num_vertices = %d;\n\n" % Nv)
  wrp.write("#declare nut_num_hole_sides = %d;\n\n" % Nr)
  write_vertices_POV(wrp, Vpos, Vlab)
  write_edges_POV(wrp, Eorg, Edst, Nr, Vpos, Vlab)
  write_faces_POV(wrp, Fnrm, Fpos, Nr)
  wrp.close()
  
  return 0

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

main()