# Implementation of module {example_path}
# Last edited on 2021-02-19 14:33:31 by jstolfi

import move
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

# MISCELLANEOUS EXAMPLE PATHS

def simple(parms):

  wd = parms['solid_raster_width']

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

  p04 = (2.75,2.5)
  p05 = (3.25,0.5)

  p06 = (4.00, 1+0*wd)
  p07 = (5.50, 1+0*wd)
  p08 = (5.50, 1+1*wd)
  p09 = (4.75, 1+1*wd)
  p10 = (4.25, 1+2*wd)
  p11 = (5.50, 1+2*wd)
  p12 = (5.50, 1+3*wd)
  p13 = (4.00, 1+3*wd)
  
  ptss = ((p01, p02, p03), (p04,p05), (p06, p07, p08, p09, p10, p11, p12, p13))
  ph = path.from_points(ptss, wd, parms)
  return ph

def scanline_fill(alt, wd, parms):

  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, wd, parms)
  mv0, dr0 = move.unpack(path.elem(ph, ix0))
  mv1, dr1 = move.unpack(path.elem(ph, ix1))
  ct = contact.make(q0, q1, mv0, mv1, parms)
  # sys.stderr.write("ph = %s ct =%s\n" % (str(ph), str(ct)))
  return ph, ct

def onion(ctr, Rc, wdc, Rf, wdf, phase, nt, parms):
  
  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, wdc, nt, phase, parms)
  phs = [ phc, ]
  Rin = Rc*cos(astep/2) - wdc/2  # Bullseye's inradius.
  Rot = 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 {Rot}.
      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
        Rot = Rot_next
      else:
        assert False
      sys.stderr.write("adding R = %.3f nt = %d\n" % (R_next,nt_next))
      phk = circle(ctr, R_next, wdf, nt_next, phase_next, parms)
      phs.append(phk)
      R_prev = R_next
      wd_prev = wdf
      phase_prev = phase_next 
      nt_prev = nt_next
  use_jumps = False
  ph = path.concat(phs, use_jumps, parms)
  return ph, Rin, Rot

def gearloose(R, zigzag, parms):

  wdfill = parms['solid_raster_width']
  wdcont = parms['contour_trace_width']
  
  Rot = R        # Outer radius of gear.
  Rin = 0.80*Rot # Inner radius of gear.
  ctr = (Rot+1,Rot+1)

  ph0 = path.from_points((ctr, ctr), wdfill, parms)
  
  ph1 = gear(ctr, Rin, Rot, 11, True, pi/7, parms)
  
  Sin = 0.20*R  # Inner radius of annulus.
  ph2 = circle(ctr, Sin, wdcont, 33, pi/11, parms)

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

  ph = path.concat((ph0, ph1, ph2, ph3, ph4), True, parms)
  return ph

def raster_rectangle(n, dp, xsz, ystep, wd, use_jumps, parms):
  # Create a list {RS} of {n} raster lines, from bottom to top:
  # Bottom-most raster oriented lef to right:
  RS = [] 
  for k in range(n):
    p = rn.add(dp, (0, k*abs(ystep)))
    q = rn.add(p, (xsz, 0))
    ph = path.from_points((p, q,), wd, parms)
    if k % 2 == 1: ph = path.rev(ph)
    RS.append(ph)
  # Sort the paths in the specified order:
  if ystep > 0:
    TS = [ RS[i] for i in range(n) ]
  else:
    TS = [ RS[n-1-i] for i in range(n) ]
  # Concatenates the paths in the proper order, with links if possble:
  P = path.concat(TS, use_jumps, parms)
  return P
  # ----------------------------------------------------------------------

# CONTOUR COMPONENTS

def circle(ctr, R, wd, nt, phase, parms):

  # 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, wd, parms)
  return ph

def gear(ctr, Rin, Rot, nt, split, phase, parms):

  wdcont = parms['contour_trace_width']

  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 = Rot
    else:
      R0 = Rot; 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, wdcont, parms)
  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, Rot, p, q, side):
  # Clips the segment {p--q} to the ring with center at the origin,
  # inner radius {Rin} and outer radius {Rot}. 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, Rot)
  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, Rot, y, side, parms):
  
  wdfill = parms['solid_raster_width']

  # Adjust the radii accounting for the round caps of the raster lines:
  Sin = Rin if Rin < 0 else Rin + wdfill/2
  Sot = Rot - wdfill/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*Rot    # Sufficient {X} distance to be outside the ring.
  p = (-Rot, y)  # On the {Y}-axis.
  q = (+Rot, 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, wdfill, parms)
 
def ring_append_clipped_seg(pts, Rin, Rot, 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 {Rot}, 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, Rot, 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, Rot):
  # Given a list {pts} of vertices of a polyline,
  # returns its intersection with the ring with center at
  # origin, inner radius {Rin}, outer radius {Rot}.
  #
  # 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, Rot, 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, Rot, y, step, phase, parms):

  wdfill = parms['solid_raster_width']
  assert step > wdfill

  # Adjust radii to account for trace round caps:
  Sin = Rin if Rin < 0 else Rin + wdfill/2
  Sot = Rot - wdfill/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 = wdfill # Height of chip removed from tooth tips.
 
  pts_unc = zigzag_segs(Rot, y, height, ytrim, phase)
  pts_clp = ring_clip_polyline(pts_unc, Sin, Sot)
  if len(pts_clp) == 0:
    return None
  else:
    return path.from_points(pts_clp, wdfill, parms)
     
def ring_fill(Rin, Rot, dmin, dmax, step, zigzag, side, parms):

  wdfill = parms['solid_raster_width']
  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, Rot, y, side, parms)
      else:
        phase = ((k + 3) % 4)/4
        assert phase == 0.0 or phase == 0.5
        phk = ring_zigzag(Rin, Rot, y, step, phase, parms)
      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")
    return None
  else:
    return path.concat(ophs, False, parms)