# Implementation of module {path_example}
# Last edited on 2021-03-21 10:29:04 by jstolfi

import path_example
import move
import move_example
import move_parms
import contact
import path
import rootray
import rootray_cart
import hacks
import rn
import sys
from math import sqrt, hypot, sin, cos, floor, ceil, inf, nan, pi

def raster_rectangle(plo, axis, n, alt, sz, step, mp_trace,mp_jump):
  
  # Create all needed rasters, links, and jumps:
  TRS = move_example.rectangle_rasters(plo, axis, n, sz, step, mp_trace)
  if alt:
    LKS0, LKS1 = move_example.rectangle_links(plo, axis, n, sz, step, mp_trace)
  else:
    JMS0, JMS1 = move_example.rectangle_jumps(plo, axis, n, sz, step, mp_jump)
  
  # Select the moves for each path:
  MVS0 = []
  MVS1 = []
  for k in range(n):
    # Append the raster line in the proper orientation:
    trk = TRS[k]
    if alt:
      # Alternating directions:
      MVS0.append(move.spin(trk, k%2))
      MVS1.append(move.spin(trk, (k+1)%2))
    else:
      # Uniform directions:
      MVS0.append(trk)
      MVS1.append(move.rev(trk))
    if k <  n-1:
      # Append the proper link or jump:
      if alt:
        # Join with links:
        if k%2 == 0:
          MVS0.append(LKS1[k])
          MVS1.append(LKS0[k])
        else:
          MVS0.append(LKS0[k])
          MVS1.append(LKS1[k])
      else:
        # Join with jumps:
        MVS0.append(JMS0[k])
        MVS1.append(JMS1[k])
        
  ph0 = path.from_moves(MVS0); path.set_name(ph0, "P0")
  ph1 = path.from_moves(MVS1); path.set_name(ph1, "P1")
  PH = [ ph0, ph1, ]
  return PH, TRS
  # ----------------------------------------------------------------------

def spiral_rectangle(pini, szx, szy, axis, mp_trace):
  wd = move_parms.width(mp_trace)
  MVS = [] # The moves
  
  p = pini          # Current point.
  ax = axis         # Axis of next move.
  d = [ szx, szy ]  # Signed length of next moves along each axis.
  k = 0 # Number of traces made.
  # Decide the displacement vectors for the next moves in each direction:
  while True:
    # Is there space for another {ax} trace:
    if abs(d[1-ax]) < 0.45*wd: break
    # Take the next move along the axis {ax} with signed length {d[ax]}:
    q = ( p[0] + d[0], p[1] ) if ax == 0 else ( p[0], p[1] + d[1] )
    mv = move.make(p, q, mp_trace)
    xd = "a" if d[ax] > 0 else "b"
    move.set_name(mv, "HV"[ax] + xd + str(k//4))
    MVS.append(mv)
    k += 1
    if k > 2:
      # Subtract the trace from the area left to cover:
      d1 = d[1-ax]
      d[1-ax] = d1 - wd if d1 > 0 else d1 + wd
    # Reverse the direction of the next {ax} trace:
    d[ax] = -d[ax]
    # Prepare for the next iteration:
    ax = 1-ax
    p = q
  ph = path.from_moves(MVS)
  pini = path.pini(ph)
  return ph
  # ----------------------------------------------------------------------

def onion(ctr, Rc, mp_cont, Rf, mp_fill, phase, nt, mp_jump):

  wdc = move_parms.width(mp_cont)
  wdf = move_parms.width(mp_fill)
  
  def reduce_nt(nt1, R1, wd1):
    # Given the number of sides {nt1} in the previous outer element (filler
    # or contour), and the radius {R1} of the vertices of the next element,
    # returns a reduced number of sides {nt1/d} for some small
    # {d}, provided that the gaps between the polygon and the circle
    # are not too large compared to the trace width {wd1}.   If it fails,
    # returns {nd1} unchanged.
    while True:
      nt1_old = nt1
      # Find a small divisor {idiv} of {nt_next}
      for idiv in 7,5,3,2:
        nt_try = nt1//idiv
        if nt1 % idiv == 0 and nt_try >= 3:
          astep_try = 2*pi/nt_try
          # Check the rel deviation {sagg} of trace edges from circle would be acceptable:
          sagg = R1*(1 - cos(astep_try/2))/wd1
          # sys.stderr.write("nt1 = %d idiv = %d sagg/wd1= %.4f\n" % (nt1, idiv, sagg/wd1))
          if sagg <= 0.20: 
            # OK, reduce {nt1} and keep trying:
            nt1 = nt_try
            sys.stderr.write("R = %.3f sagg = %.3f nt = %d\n" % (R1,sagg,nt1))
      if nt1 == nt1_old: 
        # Tried all divisors but failed to reduce:
        return nt1
    assert False
    # ----------------------------------------------------------------------

  assert Rc > 0 and wdc > 0
  assert wdf > 0
  assert nt >= 3 
  astep = 2*pi/nt
  # Start with the contour:
  phc = circle(ctr, Rc, nt, phase, mp_cont)
  phs = [ phc, ]
  Rin = Rc*cos(astep/2) - wdc/2  # Bullseye's inradius.
  Rex = Rc + wdc/2               # Bullseye's exradius.
  R_prev = Rc        # Endpoint radius of previous element.
  wd_prev = wdc      # Width of previous element
  phase_prev = phase # Phase angle of previous element.
  nt_prev = nt       # Number of sides in previous element.
  if Rf != Rc:
    while True:
      # Compute elements {R_next,nt_next,phase_next} of next circle.
      # Also update {Rin} or {Rex}.
      astep_prev = 2*pi/nt_prev # Angular increment of previous element.
      nt_next = nt_prev # May be reduced.
      ashift_prev = - ceil(2*wdf/R_prev/astep_prev)*astep_prev
      phase_next = phase_prev + ashift_prev
      if Rf < Rc:
        # Filling inwards. The radius and phase shift do not depend on {nt_next}:
        R_next = R_prev - (wd_prev + wdf)/2/cos(astep/2) # Endpoint radius of next element.
        # Try to reduce the number of segments{nt_prev}:
        nt_next = reduce_nt(nt_next, R_next, wdf)
        # Recompute the step between vertices and the phase shift:
        astep_next = 2*pi/nt_next
        # Compute inradius{Rin_next} of next element:
        Rin_next = R_next*cos(astep_next/2) - wdf/2 
        if R_next < 0 or Rin_next < Rf:
          # No space for the next element
          break
        Rin = Rin_next
      elif Rf > Rc:
        astep_next = 2*pi/nt_next
        R_next = R_prev + (wd_prev + wdf)/2/cos(astep_next/2) # Endpoint radius of next element.
        Rot_next = R_next + wdf/2  # Exradius of next element.
        if Rot_next > Rf:
          # No space for the next element
          break
        Rex = Rot_next
      else:
        assert False
      sys.stderr.write("adding R = %.3f nt = %d\n" % (R_next,nt_next))
      phk = circle(ctr, R_next, nt_next, phase_next, mp_fill)
      phs.append(phk)
      R_prev = R_next
      wd_prev = wdf
      phase_prev = phase_next 
      nt_prev = nt_next
  use_jumps = True # Don't use links between the layers.
  use_links = False
  ph = path.concat(phs, use_jumps, use_links, mp_jump)
  return ph, Rin, Rex
  # ----------------------------------------------------------------------

def gearloose(R, zigzag, mp_cont, mp_fill, mp_jump):

  wdf = move_parms.width(mp_fill)
  wdc = move_parms.width(mp_cont)
  
  Rex = R        # Outer radius of gear.
  Rin = 0.80*Rex # Inner radius of gear.
  ctr = (Rex+1,Rex+1)

  ph0 = path.from_points((ctr, ctr), mp_fill, mp_jump)
  
  ph1 = gear(ctr, Rin, Rex, 11, True, pi/7, mp_cont, mp_jump)
  
  Sin = 0.20*R  # Inner radius of annulus.
  ph2 = circle(ctr, Sin, 33, pi/11, mp_cont)

  Sot = 0.60*R  # Outer radius of annulus.
  ph3 = circle(ctr, Sot, 44, -pi/11, mp_cont)
      
  ang = 3*pi/11
  Fin =  -1 if Sin < 0 else Sin + wdc/2
  Fot = Sot - wdc/2
  dmin = -Sot
  dmax = +Sot
  step = 2*wdf if zigzag else wdf
  
  side = 0
  ph4a = ring_fill(Fin, Fot, dmin, dmax, step, zigzag, side, mp_fill, mp_jump)
  assert ph4a != None
  ph4 = path.displace(ph4a, ang, ctr)  

  use_jumps = True # Don't use links between components.
  use_links = False
  ph = path.concat((ph0, ph1, ph2, ph3, ph4), use_jumps, use_links, mp_jump)
  return ph
  # ----------------------------------------------------------------------

# MISCELLANEOUS EXAMPLE PATHS

def misc_A(mp_trace1, mp_trace2, mp_jump):

  p0 = (2,2)
  p1 = (3,3)
  p2 = (2,4)
  p3 = (1,3)

  mv01 = move.make(p0, p1, mp_trace1); move.set_name(mv01, "T0")
  mv12 = move.make(p1, p2, mp_jump)  ; move.set_name(mv12, "J0")
  mv23 = move.make(p2, p3, mp_trace2); move.set_name(mv23, "T1")

  ph = path.from_moves((mv01, mv12, mv23))
  
  return ph
  # ----------------------------------------------------------------------

def misc_B(mp_trace, mp_jump):

  p01 = (1,1)
  p02 = (1,4)
  p03 = (2,4)

  p04 = (3,5)
  p05 = (4,0)

  p06 = (5,1)
  p07 = (7,1)
  p08 = (7,2)
  p09 = (6,2)
  p10 = (5,3)
  p11 = (7,3)
  p12 = (7,4)
  p13 = (5,4)
  
  ptss = ((p01, p02, p03), (p04, p05), (p06, p07, p08, p09, p10, p11, p12, p13))
  ph = path.from_points(ptss, mp_trace, mp_jump)

  return ph
  # ----------------------------------------------------------------------

def misc_C(mp_trace, mp_jump):

  p0a = (1, 2)
  p1a = (2, 1)
  p2a = (3, 1)
  p3a = (4, 2)

  p0b = (4, 4)
  p1b = (3, 3)

  p0c = (2, 4)

  p0d = (1, 3)
  p1d = p0d

  p0e = p0c

  pha = path.from_points( ((p0a, p1a,), (p2a, p3a,)), mp_trace, mp_jump)
  phb = path.from_points( ((p0b, p1b,),), mp_trace, mp_jump)            
  phc = path.from_points( (p0c,), mp_trace, mp_jump)                    
  phd = path.from_points( ((p0d, p1d,),), mp_trace, mp_jump)            
  phe = path.from_points( ((p0e,),), mp_trace, mp_jump)                 
  
  PHS = [ pha, phb, phc, phd, phe, ]

  return PHS
  # ----------------------------------------------------------------------

def misc_D(mp_trace, mp_jump):

  TRS, JMS = move_example.misc_B(mp_trace, mp_jump)
  
  pha = path.from_moves((TRS[0], TRS[1], TRS[2],)) # Bottom is LR
  phb = path.from_moves((move.rev(TRS[0]), move.rev(TRS[3]), move.rev(TRS[2]),)) # Bottom is RL
  phc = path.from_moves((TRS[0], TRS[1], TRS[2], TRS[3],)) # Closed.

  phd = path.from_moves((TRS[4],JMS[0],TRS[5],))
  phe = path.from_moves((TRS[5],JMS[1],TRS[4],))
  
  phf = path.from_moves((TRS[6],TRS[7]))
  
  path.set_name(pha, "Pa")
  path.set_name(phb, "Pb")
  path.set_name(phc, "Pc")
  path.set_name(phd, "Pd")
  path.set_name(phe, "Pe")
  path.set_name(phf, "Pf")
  
  OPHS = [ pha, phb, path.rev(phc), phd, phe, path.rev(phf), ]
  return OPHS, TRS, JMS
  # ----------------------------------------------------------------------

def misc_E(mp_trace, mp_jump):
  
  TRS, JMS = move_example.misc_C(mp_trace, mp_jump)
  
  # Join the moves in paths:
  
  opha = path.from_moves((TRS[0],JMS[0],TRS[6],))
  ophb = path.from_moves((move.rev(TRS[2]),JMS[1],TRS[5],JMS[2],TRS[4],))
  ophc = path.rev(path.from_moves((TRS[1],JMS[3],TRS[3])))
  ophd = path.from_moves((TRS[7],JMS[4],TRS[8],))
  ophe = path.from_moves((TRS[9],))
  
  path.set_name(opha, "Pa")
  path.set_name(ophb, "Pb")
  path.set_name(ophc, "Pc")
  path.set_name(ophd, "Pd")
  path.set_name(ophe, "Pe")
  
  OPHS = [ opha, ophb, ophc, ophd, ophe, ]

  return OPHS, TRS, JMS
  # ----------------------------------------------------------------------

def misc_F(alt, mp_trace, mp_jump):

  wd = move_parms.width(mp_trace)

  p11 = (1, 1 + 0*wd) # Starting nozzle position.

  p21 = (4, 1 + 1*wd)
  p22 = (7, 1 + 1*wd)

  p31 = (4, 1 + 2*wd)
  p32 = (8, 1 + 2*wd)

  p41 = (1, 1 + 3*wd)
  p42 = (2, 1 + 3*wd)
  p43 = (4, 1 + 3*wd)
  p44 = (5, 1 + 3*wd)
  p45 = (7, 1 + 3*wd)
  p46 = (9, 1 + 3*wd)

  p51 = (1, 1 + 4*wd)
  p52 = (2, 1 + 4*wd)
  p53 = (4, 1 + 4*wd)
  p54 = (5, 1 + 4*wd)
  p55 = (7, 1 + 4*wd)
  p56 = (8, 1 + 4*wd)

  p61 = (1, 1 + 5*wd)
  p62 = (5, 1 + 5*wd)
  
  if alt:
    q0 = (1.0, 1 + 3.5*wd); ix0 =  5
    q1 = (2.0, 1 + 3.5*wd); ix1 = 15
    pts= \
      ( ( p11, ),
        ( p21, p22, p32, p31, ),
        ( p41, p42, ), ( p43, p44, ), ( p45, p46, p56, p55, ), ( p54, p53, ), ( p52, p51, p61, p62, ),
        # ( p41, p42, ), ( p43, p44, ), ( p45, p46, ), ( p54, p53, ), ( p52, p51, p61, p62, ),
      )
  else:
    q0 = (1.0, 1 + 3.5*wd); ix0 =  3
    q1 = (2.0, 1 + 3.5*wd); ix1 =  9
    pts = \
       ( ( p11, ),
        ( p21, p22, ), 
        ( p31, p32, ),
        ( p41, p42, ), ( p43, p44, ), ( p45, p46, ), 
        ( p51, p52, ), ( p53, p54, ), 
        ( p55, p56, ),
        # ( p55, p56, ),
        ( p61, p62, ),
      )
  ph = path.from_points(pts, mp_trace, mp_jump)
  mv0, dr0 = move.unpack(path.elem(ph, ix0))
  mv1, dr1 = move.unpack(path.elem(ph, ix1))
  ct = contact.make(q0, q1, mv0, mv1)
  # sys.stderr.write("ph = %s ct =%s\n" % (str(ph), str(ct)))
  return ph, ct

# CONTOUR COMPONENTS

def circle(ctr, R, nt, phase, mp_trace):

  wd = move_parms.width(mp_trace)
  mp_jump = None # Should be just traces.

  # Leave a small gap at the end for visual clarity.
  # Also it may avoid a bump on the surface.
  pts = []
  astep = 2*pi/nt
  agap = wd/2/R # Angular size of gap between the endpoints.
  aini = phase + agap/2        # Start extruding at this angle.
  afin = phase + 2*pi - agap/2 # Stop extruding at this angle.
  eps = 0.01*astep
  aant = None # Previous value of angle {a}.
  pant = None # Corresponding vertex of unclipped polygon.
  for k in range(nt+1):
    a = phase + k*astep
    p = rn.add(ctr, (R*cos(a), R*sin(a)))
    if k == 0:
      # Starts inside the gap.  Omit:
      pass
    elif aant < aini - eps and a > aini + eps:
      # Trace crosses out of the gap.
      # Add only the final part of the trace:
      r = (aini - aant)/astep # Not accurate, but will do.
      q = rn.mix(1-r, pant, r, p)
      pts.append(q)
      pts.append(p)
    elif aant < afin - eps and a > afin + eps:
      # Trace crosses into the gap.
      # Add only the initial part of the trace:
      pts.append(pant)
      r = (afin - aant)/astep # Not accurate, but will do.
      q = rn.mix(1-r, pant, r, p)
      pts.append(q)
    elif a >= aini and a <= afin:
      # Add vertexto path:
      pts.append(p)
    else:
      # Omit vertex and previous side, if any:
      pass
    pant = p
    aant = a
  ph = path.from_points(pts, mp_trace, mp_jump)
  return ph

def gear(ctr, Rin, Rex, nt, split, phase, mp_trace, mp_jump):

  wd = move_parms.width(mp_trace)

  ptss = [] # List of lists of points, each being a separate subpath.
  pts = None # Points of current subpath, or {None} between subpaths.
  # Incomplete gear:
  for k in range(2*nt):
    # Generate the half-tooth number {k}:
    if k % 2 == 0:
      R0 = Rin; R1 = Rex
    else:
      R0 = Rex; R1 = Rin
    # Compute points of gear half-tooth:
    a = (k+0.00)*pi/nt + phase
    p1 = rn.add(ctr, (R0*cos(a), R0*sin(a)))
    b = (k+0.50)*pi/nt + phase
    p2 = rn.add(ctr, (R0*cos(b), R0*sin(b)))
    p3 = rn.add(ctr, (R1*cos(b), R1*sin(b)))
    c = (k+1.00)*pi/nt + phase
    p4 = rn.add(ctr, (R1*cos(c), R1*sin(c)))
    if split and (k == 0 or k == 1 or k == 4 or k == 5):
      # Omit the half-tooth, ending the current path if any::
      if pts != None:
        ptss.append(pts)
        pts = None
    else:
      # Add the half-tooth:
      if pts == None:
        # First point of subpath:
        pts = [ p1 ]
      else:
        # Must be joined:
        assert p1 == pts[-1]
      pts.append(p2)
      pts.append(p3)
      pts.append(p4)
  if pts != None: ptss.append(pts)
  ph = path.from_points(ptss, mp_trace, mp_jump)
  return ph

# FILLING COMPONENTS

def circle_rootray(p, q, R):
  # Returns a root-ray for the straight line path {r(t) = p + t (q-p)}
  # relative to a circle with center at the origin and radius {R}.
  
  if R <= 0: return rootray_vacuum()
  # Corresponding points for the unit ball:
  pu = rn.scale(1/R, p)
  qu = rn.scale(1/R, q)
  K = rootray_cart.ball(pu,qu)
  # sys.stderr.write("R = %10.6f y = %10.6f H = %10.6f K = %s\n" % (R, y, H, str(K)))
  return K
  
def ring_clip_segment(Rin, Rex, p, q, side):
  # Clips the segment {p--q} to the ring with center at the origin,
  # inner radius {Rin} and outer radius {Rex}. The result is a tuple with
  # zero, one, or two pairs of points. Each pair is the endpoints of one
  # part of the clipped segment.
  #
  # When this procedure is used for filling, the caller must adjust the 
  # radii to account for the thickness of contours and filling traces.
  # sys.stderr.write("\n")
  Kot = circle_rootray(p, q, Rex)
  Kin = circle_rootray(p, q, Rin)
  Ksg = (+1, (0, 1,),)
  K = rootray.intersection(rootray.difference(Kot, Kin), Ksg)
  if K == None: 
    # Scan line does not meet the ring:
    return []
  assert type(K) is list or type(K) is tuple and len(K) == 2
  sK,tK = K # Starting status and root list.
  # sys.stderr.write("tK = %s\n" % str(tK))
  assert sK == +1
  nt = len(tK)
  # Decide which pair(s) of roots to use. The low indices are {ka,kb}.. 
  if nt == 0:
    # Segment is entirely outside the ring.
    return []
  elif nt == 2:
    # Segment crosses the outer circle but not the inner one:
    ka = 0; kb = None
  elif nt == 4:
    # Segment crosses both circles:
    if side < 0:
      ka = 0; kb = None
    elif side > 0:
      ka = 2; kb = None
    else: 
      ka = 0; kb = 2
  else:
    assert False, "inconsistent number of roots"
  assert ka != None
  
  # Now generate the point pairs list:
  segs = [] # List of pairs of points.
  for k in ka, kb:
    if k != None:
      ptpair = []
      for i in k,k+1:
        t = tK[i]
        pt = rn.mix(1-t, p, t, q)
        ptpair.append(pt)
      segs.append(ptpair)
  return segs
      
def ring_raster(Rin, Rex, y, side, mp_fill, mp_jump):
  
  wdf = move_parms.width(mp_fill)

  # Adjust the radii accounting for the round caps of the raster lines:
  Sin = Rin if Rin < 0 else Rin + wdf/2
  Sot = Rex - wdf/2

  # Quick tests for emptiness:
  if Sot <= 0: return None
  if Sin > 0 and Sot <= Sin: return None 
  if abs(y) >= Sot: return None

  # Choose {p,q} that define the ray:
  h = 1.1*Rex    # Sufficient {X} distance to be outside the ring.
  p = (-Rex, y)  # On the {Y}-axis.
  q = (+Rex, y)  # On the same scan-line, far away from the axis.
  
  # Find the ray intersection with the ring, 
  segs = ring_clip_segment(Sin, Sot, p, q, side)
  if len(segs) == 0: return None
  return path.from_points(segs, mp_fill, mp_jump)
 
def ring_append_clipped_seg(pts, Rin, Rex, p, q):
  # The parameter {pts} must be a list of lists of points in the
  # format suitable for {path.from_points}.
  # Clips the segment {p--q} to the ring  with center at the origin,
  # inner radius {Rin} and outer radius {Rex}, and appends 
  # the pieces to {pts}, with breaks at the gaps.
  assert type(pts) is list
  # Get last chain {cant} of {pts} and its final point {pant}:
  if len(pts) == 0:
    cant = None; pant = None
  else:
    cant = pts[-1] # Last CRS in {pts}.
    assert len(cant) >= 2
    pant = cant[-1]
  # Clip {p--q} to the ring:
  # sys.stderr.write("p = %s q = %s\n" % (str(p), str(q)))
  segs = ring_clip_segment(Rin, Rex, p, q, 0)
  # sys.stderr.write("segs = %s\n" % str(segs))
  assert type(segs) is list or type(segs) is tuple
  # Now {segs} should be a list of pairs of points.
  for seg in segs:
    # Get first point of next segment {seg}:
    assert type(seg) is list or type(seg) is tuple
    assert len(seg) == 2
    ppos = seg[0]
    if pant == None or ppos != pant:
      # Start a new CRS:
      cnew = [ seg[0], seg[1], ]
      pts.append(cnew)
      cant = cnew
    else:
      # Append point to last CRS:
      cant.append(seg[1])
    pant = seg[1]

def ring_clip_polyline(pts_unc, Rin, Rex):
  # Given a list {pts} of vertices of a polyline,
  # returns its intersection with the ring with center at
  # origin, inner radius {Rin}, outer radius {Rex}.
  #
  # The output is a list of lists of points suitable 
  # for {path;from_points}.
  assert type(pts_unc) is list or type (pts) is tuple
  # sys.stderr.write("pts_unc = %s\n" % str(pts_unc))
  pts_clp = []
  pant = None
  for p in pts_unc:
    assert hacks.is_point(p)
    if pant != None:
      ring_append_clipped_seg(pts_clp, Rin, Rex, pant, p)
    pant = p
  # sys.stderr.write("pts_clp = %s\n" % str(pts_clp))
  return pts_clp

def zigzag_segs(w, y, h, ytrim, phase):
  # Returns a list of points that are the vertices of the axes 
  # of a zigzag path with trimmed 60 degree teeth. The path 
  # is not clipped against anything.
  #
  # The medial line of the zigzag is the scanline at height {y}. 
  # The tips of the untrimmed teeth lie at ordinates {y+h} and
  # {y-h}.  The depth of trimming is {ytrim}. 
  # The polygonal line extends from at least {-w} to {+w} in {X}.
  #
  # There is the tip of a tooth at (phase*stride, y + h) where {stride} is the 
  # horizontal pitch (distance between successive tooth tips).
  assert ytrim < h 
  
  # Some geometric parameters:
  stride = 4*h/sqrt(3) # Distance between teeth on the same side of {L}.
  xtrim = ytrim/sqrt(3) # Half-width of flat part of tooth.

  # Start on {L}, on the {-X} side, well away from the circle:
  n = int(ceil(w/stride + 0.00001)) + 1 # Sufficient number of teeth.
  assert n >= 1
  phase = (phase % 1.0) # Reduce the phase to {[0 _ 1)}
  assert phase >= 0 and phase < 1.0
  xprev = (- n + phase)*stride # Ideal toot point abscissa.

  sprev = +1 # Tooth tip direction:  {+1} or {-1} for above or below {L}.
  ymin = y - h; ymax = y + h # {Y} range for points on the zigzag axis
  
  pts = []    # List of points.

  # Now generate the teeth:
  while xprev <= +w:
    # The point {(xprev,y+sprev*h} is the previous ideal tooth tip.
    # We generate the sloping part of the tooth after {xprev},
    # plus the flat part at the top of the next toot.
    
    # Compute the position {(xnext,y+snext*h)} of the next ideal tooth tip:
    xnext = xprev + stride/2
    snext = -sprev

    # Get endpoints of sloping part:
    p = (xprev + xtrim, y + sprev*(h - ytrim))
    pts.append(p)

    q = (xnext - xtrim, y + snext*(h - ytrim))
    pts.append(q)
    
    # Advance to next tooth:
    xprev = xnext
    sprev = snext
  
  assert len(pts) >= 2
  return pts

def ring_zigzag(Rin, Rex, y, step, phase, mp_fill, mp_jump):

  wdf = move_parms.width(mp_fill)
  assert step > wdf

  # Adjust radii to account for trace round caps:
  Sin = Rin if Rin < 0 else Rin + wdf/2
  Sot = Rex - wdf/2

  # Quick tests for emptiness:
  if Sot <= 0: return None
  if Sin > 0 and Sot <= Sin: return None 
  if abs(y) >= Sot: return None
  
  # Some geometric parameters:
  height = step       # Displacement of ideal tooth tips from midline.
  ytrim = wdf # Height of chip removed from tooth tips.
 
  pts_unc = zigzag_segs(Rex, y, height, ytrim, phase)
  pts_clp = ring_clip_polyline(pts_unc, Sin, Sot)
  if len(pts_clp) == 0:
    ph = None
  else:
    ph = path.from_points(pts_clp, mp_fill, mp_jump)
  return ph
     
def ring_fill(Rin, Rex, dmin, dmax, step, zigzag, side, mp_fill, mp_jump):

  if zigzag: side = 0 # Ignore {side} in zigzag fill.

  # Find range of {k}: 
  kmin = int(floor(dmin/step - 0.0001))-1 # Start below {dmin}, just to be sure.
  kmax = int(ceil(dmax/step + 0.0001))+1  # Stop above {dmax}, just to be sure
  k = kmin
  ophs = [ ] # List of paths
  while (k <= kmax):
    y = k*step
    if dmin <= y and y <= dmax:
      if not zigzag or (k % 2) == 0:
        phk = ring_raster(Rin, Rex, y, side, mp_fill, mp_jump)
      else:
        phase = ((k + 3) % 4)/4
        assert phase == 0.0 or phase == 0.5
        phk = ring_zigzag(Rin, Rex, y, step, phase, mp_fill, mp_jump)
      if phk != None:
        ophk = phk if (k % 2) == 0 else path.rev(phk)
        # sys.stderr.write("ophk =%s\n" % str(ophk))
        ophs.append(ophk)
    k = k + 1
  if len(ophs) == 0:
    sys.stderr.write("path is empty!\n")
    ph = None
  else:
    # Concatenate the elements, using links if appropriate:
    use_jumps = False
    use_links = False
    ph = path.concat(ophs, use_jumps, use_links, mp_jump)
  return ph