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 slicing_naf_vertices()\n")
  wrp.write("  #local Nv = slicing_naf_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, Vpos, Vlab):
  # Writes to {wrp} the indices of the endpoints of the edges of the fan.
  
  Nv = len(Vpos) - 1
  Ne = len(Eorg) - 1

  # Write the POV-Ray preambles:
  wrp.write("#macro slicing_naf_edges()\n")
  wrp.write("  #local Ne = slicing_naf_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 upper chain edge of plaza [0].\n")
  wrp.write("  //   200-299 inner??? chain edge of plaza [0].\n")
  wrp.write("  //   300-399 upper chain edge of plaza [1].\n")
  wrp.write("  //   400-499 inner??? chain edge of plaza [1].\n")
  wrp.write("  //   0       other.\n")
  wrp.write("  // \n")
  wrp.write("  // The index of the edge in the chain is {ty[ke] % 100}.\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):
    kv_org = Eorg[ke];
    kv_dst = Edst[ke]; 
    
    # Determine type and index of origin and destination:
    lab_org = Vlab[kv_org]; # Label of vertex {kv_org}.
    kc_org = int(re.sub(r"^v[io][.][01][.]", "", lab_org))
    
    lab_dst = Vlab[kv_dst]; # Label of vertex {kv_dst}.
    kc_dst = int(re.sub(r"^v[io][.][01][.]", "", lab_dst))
    
    vo0_pat = "^vo[.]0[.][0-9]$"      # upper vertex of plaza [0].
    vo0_org = re.match(vo0_pat, lab_org); vo0_dst = re.match(vo0_pat, lab_dst)
    
    vo1_pat = "vo[.]1[.][0-9]+$"      # upper vertex of plaza [1].
    vo1_org = re.match(vo1_pat, lab_org); vo1_dst = re.match(vo1_pat, lab_dst)
    
    vi0_pat = "^vi[.]0[.][0-9]$"      # Inner??? vertex of plaza [0].
    vi0_org = re.match(vi0_pat, lab_org); vi0_dst = re.match(vi0_pat, lab_dst)
    
    vi1_pat = "vi[.]1[.][0-9]+$"      # Inner??? vertex of plaza [1].
    vi1_org = re.match(vi1_pat, lab_org); vi1_dst = re.match(vi1_pat, lab_dst)

    if vo0_org and vo0_dst:
      # Edge is upper chain of plaza [0]:
      ty = 100 + kc_org
    elif vi0_org and vi0_dst:
      # Edge is inner??? chain of plaza [0]:
      ty = 100 + kc_org
    elif vo1_org and vo1_dst:
      # Edge is upper chain of plaza [1]:
      ty = 300 + kc_org
    elif vi1_org and vi1_dst:
      # Edge is inner??? chain of plaza [1]:
      ty = 400 + kc_org
    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, f"edge {lab_org}--{lab_dst} too short"
    
  # Write the POV-Ray postamble:
  wrp.write("  E\n")
  wrp.write("#end\n\n")
 
def write_faces_POV(wrp, Fnrm, Fpos):
  # Writes the face planes to {wrp} as POV-ray.
  
  Nf = len(Fnrm) - 1

  # Write the POV-Ray preambles:
  wrp.write("#macro slicing_naf_faces()\n")
  wrp.write("  #local Nf = slicing_naf_num_faces;\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 order is upper chain faces, inner??? chain faces.\n")
  wrp.write("  // inner???-upper connecting faces,\n")
  wrp.write("  // and finally the plazas.\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")


  Fnrm, Fpos, Eorg, Edst = write_faces_OBJ(wro, Pind, Vpos)
  wro.close()
  
  Nv = len(Vpos) - 1
  Ne = len(Eorg) - 1
  Nf = len(Fnrm) - 1
  
  wrp = open(file_pref + ".inc", 'w')
  wrp.write("#declare slicing_naf_upper_chain_even_radius = %.*f;\n\n" % (prec, Ru))
  wrp.write("#declare slicing_naf_upper_chain_odd_radius = %.*f;\n\n" % (prec, Ro_odd???))
  wrp.write("#declare slicing_naf_inner???_chain_even_radius = %.*f;\n\n" % (prec, Rieven???))
  wrp.write("#declare slicing_naf_inner???_chain_odd_radius = %.*f;\n\n" % (prec, Riodd???))
  wrp.write("#declare slicing_naf_thickness = %.*f;\n\n" % (prec, Hx))
  wrp.write("#declare slicing_naf_min_elevation??? = %.5f; // Degrees.\n\n" % Amin???)
  wrp.write("#declare slicing_naf_max_elevation??? = %.5f; // Degrees.\n\n" % Aspoke)
  
  wrp.write("#declare slicing_naf_num_faces = %d;\n\n" % Nf)
  wrp.write("#declare slicing_naf_num_edges = %d;\n\n" % Ne)
  wrp.write("#declare slicing_naf_num_vertices = %d;\n\n" % Nv)
  wrp.write("#declare slicing_naf_num_upper_chain_edges = %d;\n\n" % Ni)
  wrp.write("#declare slicing_naf_num_inner???_chain_edges = %d;\n\n" % Nei???)
  wrp.write("#declare slicing_naf_upper_chain_shape = %d;\n\n" % oshape???)
  write_vertices_POV(wrp, Vpos, Vlab)
  write_edges_POV(wrp, Eorg, Edst, Vpos, Vlab)
  write_faces_POV(wrp, Fnrm, Fpos)
  wrp.close()

    
  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)


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

  wro.write("\n")

  return Pind, Vpos, Vlab

def map_plaza(P, func):
  Ni = len(P['vo']) - 1;
  Nei??? = len(P['vi']) - 1;
  Pmap = make_null_plaza(Ni,Nei???)

  for kc in range(Ni+1):
    Pmap['vo'][kc] = func(P['vo'][kc])

  for kc in range(Nei???+1):
    Pmap['vi'][kc] = func(P['vi'][kc])
    
  return Pmap

def make_null_plaza(Ni, Nei???):
  # Creates an array like the total vertex table but with all
  # entries set to {None}.
  Pnull = {}
  Pnull['vo'] = [ None ]*(Ni+1)
  Pnull['vi'] = [ None ]*(Nei???+1)
  return Pnull


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 write_faces_OBJ(wro, Pind, Vpos):
  # Writes the faces of the fan object {to OBJ file {wro}, assuming that
  # the vertex indices are in {Pind}.
  
  # 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}.
  
  Ni = len(Pind[0]['vo']) - 1;   # Number of upper edges per plaza.
  Nei??? = len(Pind[0]['vi']) - 1;   # Number of inner??? edges per plaza.

  Nv = len(Vpos) - 1              # Expected total number of vertices.
  Ne = 3*(Ni + Nei??? + 2)          # Expected total number of edges.
  Nf = 2 + (Ni + Nei??? + 2)        # Expected total number of face planes.
  assert Nv == 2*(Ni + Nei??? + 2), "Nv,Ni,Nei??? inconsistent"
  
  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, debug = False):
    # Title for a face or set of faces.
    nonlocal debug_face, Fv, Fn
    wro.write("\n");
    wro.write("# %s\n" % tit);
    sys.stderr.write("writing %s \n" % (tit));
    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

  Ptit("upper chain faces")
  for kc in range(Ni):
    Bof()
    Prtv(Pind[0]['vo'][kc])
    Prtv(Pind[1]['vo'][kc])
    Prtv(Pind[1]['vo'][kc+1])
    Prtv(Pind[0]['vo'][kc+1])
    Eof()

  Ptit("inner??? chain faces")
  for kc in range(Nei???):
    Bof()
    Prtv(Pind[1]['vi'][kc])
    Prtv(Pind[0]['vi'][kc])
    Prtv(Pind[0]['vi'][kc+1])
    Prtv(Pind[1]['vi'][kc+1])
    Eof()
    
  # The inner???-upper connecting faces must come after the chain faces:

  Ptit("bottom inner???-upper connecting face")
  Bof()
  Prtv(Pind[0]['vi'][0])
  Prtv(Pind[1]['vi'][0])
  Prtv(Pind[1]['vo'][0])
  Prtv(Pind[0]['vo'][0])
  Eof()

  Ptit("top inner???-upper connecting face")
  Bof()
  Prtv(Pind[0]['vo'][Ni])
  Prtv(Pind[1]['vo'][Ni])
  Prtv(Pind[1]['vi'][Nei???])
  Prtv(Pind[0]['vi'][Nei???])
  Eof()
  
  # The plazas should come last:
  # Beware that if {oshape???==1} or {oshape???==2} the first two upper edges are concave
  # so must start with a spoke edge:
   
  Ptit(f"plaza[0]")
  Bof()
  Prtv(Pind[0]['vi'][0])
  for kc in range(Ni+1):
    Prtv(Pind[0]['vo'][kc])
  for kc in range(Nei???):
    Prtv(Pind[0]['vi'][Nei???-kc])
  Eof()
  
  Ptit(f"plaza[1]")
  Bof()
  Prtv(Pind[1]['vi'][Nei???])
  for kc in range(Ni+1):
    Prtv(Pind[1]['vo'][Ni-kc])
  for kc in range(Nei???):
    Prtv(Pind[1]['vi'][kc])
  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