#! /usr/bin/python3
# Last edited on 2024-06-09 10:18:08 by stolfi

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

# The program takes three parameters from the command line: an edge count {Neo},
# a general shape code {oshape}, and a boolean flag {inner}.
#
# In any case, the object is a fan-like solid that has two large
# congruent faces (the /plazas/) perpendicular to the {X}-axis,
# connected by a ring of rectangular faces. The two plazas are indexed
# by {km} in {0,1}; plaza [0] has positive {X} and plaza [1] is its
# mirror image across the plane {X = 0}.
#
# The border of each plaza consists of two poligonal chains, the /outer
# chain/ and the /inner chain/, connected by two edges (the /spokes/)
# whose supporting lines intercept the {X}-axis.
#
# The parameter {Neo} specifies the number of edges in the outer chain.
# If the {inner} parameter is false, the inner chain is reduced to a
# single vertex ({Nei == 0}). Otherwise the inner chain has the same
# number of edges as the outer one ({Nei == Neo}). The outer chain has
# two /extreme vertices/ and {Neo-1}, /interior vertices/, and similarly
# for the innner chain.
#
# A set of points is said to be /cocircular/ is they lie on the same circle
# parallel to the {YZ} plane.
#
# The /elevation/ of a vertex or spoke is its angular coordinate around
# the {X}-axis, measured CCW from {Z=0} plane as seen from {(+oo,0,0)}.
#
# Currently the allowed values of {oshape} parameter may be:
#
#   0 The number {Neo} must be a power of 2. All the {Neo+1} vertices of
#     the outer chain are convex and cocircular with center on the
#     {X}-axis. If {inner} is true, the inner chain too has {Neo+1}
#     vertices, also cocircular, with all interior vertices concave. The
#     spoke elevations {Amin,Amax} are strictly between 0 and 90
#     degrees.
#   
#   1 The number {Neo} must be a power of 2 plus 2. Of the {Neo+1}
#     vertices of the outer chain, {Neo/2} have shallow concave angles
#     and alternate with {Neo/2+1} convex ones. Each subset is
#     cocircular with center on the {X}-axis. If {inner} is true, the
#     inner chain has {Neo} edges, all cocircular with center on the
#     {X}-axis, and all its interior vertices concave. The spoke
#     elevations {Amin,Amax} are strictly between 0 and 90 degrees.
#
#   2 The number {Neo} must be even. The the extremal vertices of the
#     outer chain are convex, all its {Neo-1} interior angles are concave.
#     The vertices are cocircular, but the center is not on the
#     {X}-axis; in fact the outer chain is the same as that for {oshape=0},
#     but is mirrored across its chord line. If {inner} is true, the inner
#     chain too has {Neo} edges, with all vertices
#     cocircular with center on the {X}-axis, and the interior vertices
#     are all concave. The elevations {Amin,Amax} are symmetric about the
#     vertical ({Amax = 180 - Amin}), with {Amin} are strictly between 0
#     and 90 degrees.
#
#   3 The number {Neo} must be a power of 2 plus 2.  Of its {Neo+1} vertices,
#     {Neo/2} have deep concave angles and alternate with {Neo/2+1} 
#     convex ones. Each subset is cocircular
#     with center on the {X}-axis. If {inner} is true, the inner also
#     chain has {Neo} edges; {Neo/2} of the vertices are concave,
#     the other {Neo/2+1} are convex, with each subset cocircular with
#     center on the {X}-axis. The spoke elevations {Amin,Amax} are
#     symmetric about the vertical ({Amax = 180 - Amin}), with {Amin}
#     are strictly between 0 and 90 degrees.

def ref_plaza_vertices(Neo, Ro_even, Ro_odd, Ri_even, Ri_odd, Amin, Amax, inner, oshape):
  # Computes the vertices of a plaza at {X=0}, as an array of 
  # pairs {(Yk, Zk)}. Parameters:
  
  # {Neo}      Number of edges in the outer chain.
  # {Ro_even}  Circumradius of even vertices in outer chain in each plaza.
  # {Ro_odd}   Circumradius of odd vertices in outer chain in each plaza.
  # {Ri_even}  Circumradius of even vertices in inner chain in each plaza (0 if not {inner}).
  # {Ri_odd}   Circumradius of odd vertices in inner chain in each plaza (0 if not {inner}).
  # {Amin}     Min spoke elevation (degrees).
  # {Amax}     Min spoke elevation (degrees).
  # {Rc}       Radius of chamfer vertices on top face.
  # {Ha}       Actual distance between plazas.
  # {inner}    Bool: should the inner chain more than one vertex?
  # {oshape}   Shape of outer chain.

  sys.stderr.write("ref_plaza_vertices:\n")
  sys.stderr.write("  Neo = %d  oshape = %d  inner = %d \n" % (Neo, oshape, int(inner)))
  sys.stderr.write("  Amin = %.4f  Amax = %.4f\n" % (Amin, Amax))
  sys.stderr.write("  Ro_even = %.*f  Ro_odd = %.*f\n" % (prec, Ro_even, prec, Ro_odd))
  sys.stderr.write("  Ri_even = %.*f  Ri_odd = %.*f\n" % (prec, Ri_even, prec, Ri_odd))
  
  if inner:
    assert Ri_even > 30*eps and Ri_odd >= Ri_even, "{Ri_even,Ri_odd} should be positive"
  else:
    assert Ri_even == 0 and Ri_odd == 0, "{Ri_even,Ri_odd} should be negative"
  
  # The vertices are returned as a {dict} with keys {'vo','vi'},
  # each an array of points of {\RR^3}. Both chains are stored in CCW order
  # around the origin as seen from {(+oo,0,0)}, and indexed from 0.

  # Spoke elevations in radians:
  ang_min = radians(Amin)         # Elevation of lower spoke.
  ang_max = radians(Amax)         # Elevation of upper spoke.
  ang_mid = (ang_min + ang_max)/2 # Mean elevation of outer and inner vertices.

  # Compute outer chain vertices, initially cocircular with center on {X}-axis:
  assert Ro_odd > Ri_odd + 30*eps and Ro_even >= Ro_odd, "{Ro_even,Ro_odd} too small"
  vo = [ None ]*(Neo+1)
  dango = (ang_max - ang_min)/Neo  # Elevation increment per outer edge.
  mc = ( Ro_even*cos(ang_mid), Ro_even*sin(ang_mid) ) # Middle of chord.
  dv, ev = rn.dir(mc) # Direction and distance of middle of chord.
  ev = Ro_even*cos((ang_max - ang_min)/2) # Distance from origin to middle of chord.
  for iv in range(Neo+1):
    Ro = Ro_even if iv % 2 == 0 else Ro_odd;
    ai = ang_min + iv*dango  # Elevation of vertex.
    Yi = Ro*cos(ai)           # {Y} coord of {vo[]iv]}.
    Zi = Ro*sin(ai)           # {Z} coord of {vo[]iv]}.
    vo[iv] = ( Yi, Zi )
    if oshape == 2:
      # Mirror vertex about its chord:
      dot = rn.dot(vo[iv], dv)
      vo[iv] = rn.mix(1, vo[iv], 2*(ev-dot), dv)
  
  # Compute inner chain vertices, initially cocircular with center on {X}-axis:
  # Inner chain has {Nei} edges:
  if inner:
    Nei = Neo
  else:
    Nei = 0
  vi = [ None ]*(Nei+1)
  if inner:
    dangi = (ang_max - ang_min)/Nei  # Elevation increment per inner edge.
    for iv in range(Nei+1):
      Ri = Ri_even if iv % 2 == 0 else Ri_odd
      ai = ang_min + iv*dangi      # Elevation of vertex.
      Yi = Ri*cos(ai)              # {Y} coord of {vi[kc]}.
      Zi = Ri*sin(ai)              # {Z} coord of {vi[kc]}.
      vi[iv] = ( Yi, Zi )
  else:
    # Only one vertex on {X} axis:
    vi[0] = ( 0, 0 )
   
  Pref = { 'vo': vo, 'vi': vi }
  return Pref
  
def all_vertices(Pref, Ha):
  # The structure {Pref} must be a {dict} that gives the {Y,Z} coordinates of the vertices of the
  # reference plaza of the object, as produced by {ref_plaza_vertices}.
  # Returns an array {Ptot[0..1]} where each element is a copy
  # of {Pref} properly shifted in {X} by {±Ha/2}.

  sys.stderr.write("Ha = %.*f\n" % (prec, Ha))

  Ptot = [ None, None ]

  for km in range(2):
  
    def shift_point(p):
      nonlocal km, Ha;
      return ( (1-2*km)*Ha/2, p[0], p[1] )

    Ptot[km] = map_plaza(Pref, shift_point)
  return Ptot

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

  for kc in range(Neo+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(Neo, Nei):
  # Creates an array like the total vertex table but with all
  # entries set to {None}.
  Pnull = {}
  Pnull['vo'] = [ None ]*(Neo+1)
  Pnull['vi'] = [ None ]*(Nei+1)
  return Pnull
  
def write_vertices_OBJ(wro, Ptot):
  # Expect {Ptot} to be a an array of 2 plazas comprising the
  # vertices of a fan object. Writes the vertices to {wro} in OBJ
  # format.
  # 
  # Rounds each coordinate to even multiple of {eps}. Returns an array
  # {Pind} with the same shape as {Ptot} 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.  The vertices are listed in ccw order around the origin 
  # as seen from {(+oo,0,0)}, first all on plaza [0], then all on plaza [1].
  # 
  # Also returns a list {Vlab} where element {Vlab[i]} is the
  # position in {Ptot} of vertex {Vpos[i]}. Namely
  # "vo.0.5" to mean {Ptot[0]['vo'][5]} or "vi.1.4" to mean {Ptot[1]['vi'][4]}.

  Neo = len(Ptot[0]['vo']) - 1  # Number of outer edges .
  Nei = len(Ptot[0]['vi']) - 1  # Number of inner edges .

  Nv = 2*(Neo + Nei + 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)
    
  Pind = [ None, None ]

  kv = 0;
  for km in range(2):
    Pind[km] = make_null_plaza(Neo,Nei)
    for kc in range(Neo+1):
      kv += 1;
      Prtv(Ptot[km]['vo'][kc], f"vo.{km}.{kc}")
      Pind[km]['vo'][kc] = kv
    for kc in range(Nei+1):
      kv += 1;
      Prtv(Ptot[km]['vi'][kc], f"vi.{km}.{kc}")
      Pind[km]['vi'][kc] = kv

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

  wro.write("\n")

  return Pind, Vpos, Vlab
  
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}.
  
  Neo = len(Pind[0]['vo']) - 1;   # Number of outer 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*(Neo + Nei + 2)          # Expected total number of edges.
  Nf = 2 + (Neo + Nei + 2)        # Expected total number of face planes.
  assert Nv == 2*(Neo + Nei + 2), "Nv,Neo,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("outer chain faces")
  for kc in range(Neo):
    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-outer connecting faces must come after the chain faces:

  Ptit("bottom inner-outer 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-outer connecting face")
  Bof()
  Prtv(Pind[0]['vo'][Neo])
  Prtv(Pind[1]['vo'][Neo])
  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 outer edges are concave
  # so must start with a spoke edge:
   
  Ptit(f"plaza[0]")
  Bof()
  Prtv(Pind[0]['vi'][0])
  for kc in range(Neo+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(Neo+1):
    Prtv(Pind[1]['vo'][Neo-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
  
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 fan_vertices()\n")
  wrp.write("  #local Nv = fan_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 fan_edges()\n")
  wrp.write("  #local Ne = fan_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 chain edge of plaza [0].\n")
  wrp.write("  //   200-299 inner chain edge of plaza [0].\n")
  wrp.write("  //   300-399 outer 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]$"      # Outer 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]+$"      # Outer 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 outer 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 outer 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 fan_faces()\n")
  wrp.write("  #local Nf = fan_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 outer chain faces, inner chain faces.\n")
  wrp.write("  // inner-outer 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")

def main():
  Neo = int(sys.argv[1])          # Parameter defining chain lengths.
  oshape = int(sys.argv[2])      # Shape of outer chain.
  inner = int(sys.argv[3]) != 0  # Is there an inner chain?.
  
  sys.stderr.write("slicing_gan_example.py: ")
  sys.stderr.write(f" Neo = {Neo} oshape = {oshape} inner = {inner}\n")

  Ro_even = 100                  # Circumradius of outer chain.
  Hn = 10                        # Nominal thickness of fan.
    
  assert Neo >= 2 and Neo % 2 == 0, "{Neo} must be even and 2 or more"
    
  # Number of edges in inner chain:
  Nei = Neo if inner else 0

  # Choose the elevations of the spokes (degrees):
  if oshape == 0 or oshape == 1:
    Amin = 20
    Amax = 70
  elif oshape == 2 or oshape == 3:
    Amin = 60
    Amax = 120
  ang_min = radians(Amin)
  ang_max = radians(Amax)
  ang_step = (ang_max - ang_min)/Neo
  
  # Define {Ro_odd}:
  dango = (ang_max - ang_min)/Neo
  if oshape == 0:
    # Convex outer chain.
    Ro_odd = Ro_even
  elif oshape == 1:
    # Slightly dented upper chain:
     # Choose {Ro} so that the plazas are non-convex but still monotonic.
    oscale_min = cos(dango)
    oscale_max = sin(ang_max - 2*dango)/sin(ang_max - dango)
    oscale = (oscale_min + oscale_max)/2
    Ro_odd = Ro_even * oscale
  elif oshape == 2:
    # Concave upper chain:
    Ro_odd = Ro_even # Before flipping.
  elif oshape == 3:
    # Deeply dented upper chain:
    # Choose {Ro_odd} so that the plazas are not monotonic:
    oscale = cos(dango) - sqrt(3)*sin(dango)
    Ro_odd = Ro_even * oscale
  else:
    assert False, "invalid {oshape}"
  
  # Choose the inner radiii {Ri_even,Ri_odd}:
  if not inner:
    Ri_even = 0
    Ri_odd = 0
  else:
    if oshape == 0 or oshape == 1:
      # Convex or slightly dented upper chain, concave lower chain, {Amax} less than 90.
      # The highest inner vertex must be A LITTLE below the lowest outer one,
      # because the slicing planes will go through the upper chain faces.
      # The lowest outer vertex is the first one:
      Zomin = Ro_even*sin(ang_min)
      # The highest inner vertex is the last one:
      Ri_even = 0.85*Zomin/sin(ang_max)
      Ri_odd = Ri_even
    elif oshape == 2:
      # Concave upper and lower chains, {Amin,Amax} symmetrical.
      # The highest inner vertex must be WELL below the lowest outer one,
      # because the slicing planes will go between the two chains.
      # The lowest outer vertex is the middle one:
      assert Neo % 2 == 0, "{Neo} is not even"
      Zomin = Ro_odd*(1 - 2*sin(ang_step))
      # The highest inner vertex is the middle one:
      assert Nei % 2 == 0, "{Nei} is not even"
      Ri_even = 0.50*Zomin
      Ri_odd = Ri_even
    elif oshape == 3:
      # Highly dented upper and lower chains.
      # The highest inner vertex must be WELL below the lowest outer one,
      # because the slicing planes will go between the two chains.
      # The lowest outer vertex is the second one, a dent:
      Zomin = Ro_odd*sin(ang_min + ang_step)
      sys.stderr.write(f"Zomin = {Zomin:.4f}\n") 
      # The highest inner vertex is the middle one, a dent:
      assert Nei % 2 == 0, "{Nei} is not even"
      assert Nei == Neo, "{Nei,Neo} differ"
      Ri_odd = 0.50*Zomin
      Ri_even = Ri_odd*(cos(ang_step) - sqrt(3)*sin(ang_step))
    else:
      assert False, "invalid {oshape}"

  Ha = Hn - 4*eps
  
  Pref = ref_plaza_vertices(Neo, Ro_even,Ro_odd, Ri_even,Ri_odd, Amin,Amax, inner, oshape)
  assert len(Pref['vo']) - 1 == Neo
  assert len(Pref['vi']) - 1 == Nei

  Ptot = all_vertices(Pref, Ha)

  file_pref = f"out/fan_{Neo:08d}_osh{oshape:1}_ich{int(inner):1}"

  wro = open(file_pref + ".obj", 'w')
  Pind, Vpos, Vlab = write_vertices_OBJ(wro, Ptot)
  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 fan_outer_chain_even_radius = %.*f;\n\n" % (prec, Ro_even))
  wrp.write("#declare fan_outer_chain_odd_radius = %.*f;\n\n" % (prec, Ro_odd))
  wrp.write("#declare fan_inner_chain_even_radius = %.*f;\n\n" % (prec, Ri_even))
  wrp.write("#declare fan_inner_chain_odd_radius = %.*f;\n\n" % (prec, Ri_odd))
  wrp.write("#declare fan_thickness = %.*f;\n\n" % (prec, Ha))
  wrp.write("#declare fan_min_elevation = %.5f; // Degrees.\n\n" % Amin)
  wrp.write("#declare fan_max_elevation = %.5f; // Degrees.\n\n" % Amax)
  
  wrp.write("#declare fan_num_faces = %d;\n\n" % Nf)
  wrp.write("#declare fan_num_edges = %d;\n\n" % Ne)
  wrp.write("#declare fan_num_vertices = %d;\n\n" % Nv)
  wrp.write("#declare fan_num_outer_chain_edges = %d;\n\n" % Neo)
  wrp.write("#declare fan_num_inner_chain_edges = %d;\n\n" % Nei)
  wrp.write("#declare fan_outer_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()
  
  return 0

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

main()