# Implementation of module {raster_example}
# Last edited on 2021-06-04 04:02:00 by jstolfi

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

def patch_array(cols, rows, nx,ny, islands, mp_cont, mp_fill, mp_link):

  colper = 4 if islands else 3  # Periodicity of columns.
  
  wdf = move_parms.width(mp_fill)
  
  assert type(rows) is int and rows >= 1
  assert type(cols) is int and cols % colper == 1
  assert type(nx) is int and nx >= 1
  assert type(ny) is int and ny >= 5
  
  org = (2*wdf, 2*wdf)  # Coordinate of bottom corner of filling.
  xdir = (1, 0); ydir = (0, 1)          # Axis aligned rasters.
  
  # Dimensions in multiples of {wdf} assuming traces have width 0:
  nx_gp = 3 if islands else nx + 2 # Width of gap between roads.
  ixstep = nx + nx_gp              # Scan-column step between array columns.
  iystep = ny + 1                  # Scan-line step between array rows.
  nxtot = (cols-1)*ixstep + nx     # Total width of everything.
  nytot = rows*iystep              # Total height of everything.
  
  def make_raster(iy, row, col, irp, igr):
    # If {col} and {irp} are not {None}, creates a raster that is part
    # of the patch on row {row} and column {col} of the array;
    # specifically, the raster on patch scanline {irp} (in {0..ny-1}).
    # If {col} and {irp} are none, creates instead a beam spanning the
    # whole array, assumed to be at the bottom of row {row}. 
    #
    # In either case, the raster will be on global scanlne {iy} (in
    # {0..nytot}) and will be assigned to group {igr}. The contact lists
    # and link lists will be cleared.
    assert (col == None) == (irp == None)
    assert 0 <= row and row <= rows
    if col != None:
      ix0 = col*ixstep
      ix1 = ix0 + nx  
      name = "%d.%d.%d" % (row,col,irp)
    else:
      ix0 = 0
      ix1 = nxtot
      name = "%d" % row
    p0 = rn.mix(1, org, wdf, (ix0, iy))
    p1 = rn.mix(1, org, wdf, (ix1, iy))
    rmv = move.make(p0, p1, mp_fill)
    move.set_name(rmv, "T" + name)
    rph = path.from_moves((rmv,))
    path.set_name(rph, "P" + name)
    path_hp.set_group(rph, igr)
    path_hp.clear_contacts(rph)
    path_hp.clear_links(rph)
    return rph
    # ......................................................................
  
  OPHS = [ ]
  CTS = [] 
  rph_prev = [None]*cols # Previous raster in each column, possibly a beam.
  
  def make_bottom_beam(row):
    # Creates the beam at the bottom of array row {row} (in {0..cols}).
    # Adds the beam to {OPHS} and saves it into all entries
    # of {rph_prev}.  Adds simple link paths between the beam and the
    # patches below it.
    iy = row*iystep
    bph = make_raster(iy, row, None, None, 0)  # The beam at the bottom of the current row.
    
    # Create links between {bph} to rasters below it:
    for iend in 0, 1:
      col_below = (0, cols-1)[iend] # Coumn index of patch below end {iend}
      rph_below = rph_prev[col_below]
      if rph_below != None:
        name = "%d.%d.%d" % (iend,col_below,iy)
    # Add to list:
    OPHS.append(bph)

    # Reset {rph_prev}
    for col in range(cols):
      rph_prev[col] = bph
    return
    # ......................................................................

  def make_patches(row):
    # Creates the rasters of all patches of row {row} of the array.
    # Assumes that all entries {rph_prev[0..cols-1]]} are set 
    # to the beam raster just below that row. 
    # Adds the patch rasters to {OPHS} and updates {rph_prev}
    # to be the rasters at the top of each patch, or {None}.
    
    # Enumberate patches of this row:
    for col in range(cols):
      # Enumerate the scan-lines of this array slot:
      for irp in range(ny):
        iy = row*iystep + irp + 1   # Scanline index 
        # Decide if the patch on column {col} extends into this scanline:
        rc = col % colper
        island_patch = islands and rc == 2
        skip_bot = (island_patch or rc == colper-1) and (irp == 0 or irp == 1)
        skip_top = (rc == 1 or island_patch) and (irp == ny-2 or irp == ny-1)
        skip = skip_bot or skip_top
        if not skip:
          # Determine the group index {igr} of this patch:
          igr = len(OPHS)
          # Create the raster path with index {irp} in this patch:
          rph_this = make_raster(iy, row, col, irp, igr)
          OPHS.append(rph_this)
          # Remember this patch raster:
          rph_prev[col] = rph_this
        else:
          # Skip raster:
          rph_prev[col] = None
    return
    # ----------------------------------------------------------------------
   
  # Create the rasters:
  for row in range(rows):
    # Beam at bottom of row:
    make_bottom_beam(row)
    make_patches(row)
  make_bottom_beam(rows)
  
  # Add links:
  LKS = raster.create_all_raster_raster_links(OPHS, xdir, ydir, mp_link)
  
  # Add contacts:
  CTS = raster.create_all_raster_raster_contacts(OPHS, xdir, ydir)

  # Create the contours:
  OCRS = []
  
  def contour_point(ix, iy, xoff, yoff):
    # Creates a point at scan-grid point {(ix,xy)} wth offset {(xoff,yoff)}
    return rn.mix3(1, org, wdf, (ix,iy), 1, (xoff,yoff))
    # ......................................................................
  
  # Outer contour:
  wdc = move_parms.width(mp_cont)
  off = (wdc + wdf)/2  # Offset of contour from filling. 
  p0 = contour_point(    0,     0, -off, -off)
  p1 = contour_point(nxtot,     0, +off, -off)
  p2 = contour_point(nxtot, nytot, +off, +off)
  p3 = contour_point(    0, nytot, -off, +off)
  ocro = path.from_points((p0, p1, p2, p3, p0,), mp_cont, None)
  path.set_name(ocro, "CRo")
  OCRS.append(ocro)
  
  # Contours of holes and patches:
  for row in range(rows):
    for cg in range(cols//colper):
      # Contours of holes and islands in row {row}, between patches in columns {4*cg} and {4*cg+4}:
      iy0 = row*iystep
      iy3 = iy0 + ny + 1
      iy1 = iy0 + 3
      iy2 = iy3 - 3
      
      ix0 = colper*cg*ixstep + nx
      ix1 = ix0 + nx_gp
      ix2 = ix0 + ixstep
      ix6 = ix0 + (colper-1)*ixstep
      ix7 = ix1 + (colper-1)*ixstep
      ix5 = ix7 - ixstep
      
      # Hole contour:
      h00 = contour_point(ix0, iy0, +off, +off)
      h01 = contour_point(ix1, iy0, -off, +off)
      h02 = contour_point(ix1, iy2, -off, +off)
      h03 = contour_point(ix2, iy2, +off, +off)
      h04 = contour_point(ix2, iy0, +off, +off)
      h05 = contour_point(ix7, iy0, -off, +off)
      h06 = contour_point(ix7, iy3, -off, -off)
      h07 = contour_point(ix6, iy3, +off, -off)
      h08 = contour_point(ix6, iy1, +off, -off)
      h09 = contour_point(ix5, iy1, -off, -off)
      h10 = contour_point(ix5, iy3, -off, -off)
      h11 = contour_point(ix0, iy3, +off, -off)
      ocrh = path.from_points((h00,h01,h02,h03,h04,h05,h06,h07,h08,h09,h10,h11,h00,), mp_cont, None)
      path.set_name(ocrh, "CRh.%d.%d" % (row,4*cg))
      OCRS.append(ocrh)
      
      if islands:
        # Island contour:
        ix3 = ix1 + ixstep
        ix4 = ix2 + ixstep

        d0 = contour_point(ix3, iy1, -off, -off)
        d1 = contour_point(ix4, iy1, +off, -off)
        d2 = contour_point(ix4, iy2, +off, +off)
        d3 = contour_point(ix3, iy2, -off, +off)

        ocrq = path.from_points((d0,d1,d2,d3,d0,), mp_cont, None)
        path.set_name(ocrq, "CRd.%d.%d" % (row,4*cg+2))
        OCRS.append(ocrq)
  
    path.compute_contour_nesting(OCRS)

  return OCRS, OPHS, CTS
  # ----------------------------------------------------------------------

# SAMPLE PATHS FOR THE PAPER

def raster_figure_contour(mp, xLo,yLo, xS,yS,RS, yHi, xB0,xB1,yB,RB, xD,yD, off):
  # Returns the list {PCS} of contour paths for {raster_figure}.  See the comments in that
  # function for the meaning of the parameters.
  #
  PCS = [] # List of contour paths.

  MVS = []  # Moves of the current path.
  p = None  # Current point of the current contour path.
  
  tag = "c"
  p = ( xD - off, yLo - off )
  p = path.move_to(MVS, p, ( xS, yLo - off ), mp, tag)
  na = 12 # Number of traces in semicircle.
  for ia in range(na):
    a = pi*((ia + 1.0)/na - 0.5)
    q = (xS + (RS+off)*cos(a), yS + (RS+off)*sin(a))
    p = path.move_to(MVS, p, q, mp, tag)
  
  p = path.move_to(MVS, p, ( xLo - off, yHi + off ), mp, tag)
  p = path.move_to(MVS, p, ( xLo - off, yD - off ), mp, tag)
  p = path.move_to(MVS, p, ( xD - off, yD - off ), mp, tag)
  p = path.move_to(MVS, p, ( xD - off, yLo - off ), mp, tag)

  path.finish(PCS, MVS, tag)
  
  # Build the contours of the two holes:
  for ih in range(2):
    tag = "b%d" % ih
    MVS = []  # Traces of the current contour path.
    p = None       # Current point of the current contour path.
    xB = (xB0, xB1)[ih]  # {X} coord of hole center.
    nb = 8        # Number of traces in the hole contour.
    for ib in range(nb + 1):
      b = 2*pi*(ib/nb)
      q = (xB + (RB-off)*cos(b), yB + (RB-off)*sin(b))
      p = path.move_to(MVS, p, q, mp, tag)
    path.finish(PCS, MVS, tag)
    
  return PCS
  # ----------------------------------------------------------------------

def raster_figure_filling(mp_fill, mp_jump, xLo,yLo, xS,yS,RS, yHi, xB0,xB1,yB,RB, xD,yD, wd, version):
  # Returns the lists {PFS} and {PLS} for {raster_figure}.  See the comments in that
  # function for the meaning of the parameters.
  #

  # Define the endpoints of the rasters:
  org = (xD - 2*wd, yLo - 2*wd)  # Starting nozzle position.

  y0 = yLo + 0*wd
  dxS0 = 0
  p00 = (xD, y0)
  p01 = (xS + dxS0, y0)

  y1 = yLo + 1*wd
  dxS1 = sqrt(max(0, RS**2 - (yS-y1)**2))
  p10 = (xD, y1)
  p11 = (xS + dxS1, y1)

  y2 = yLo + 2*wd
  dxB2 = sqrt(max(0, RB**2 - (yB-y2)**2))
  dxS2 = sqrt(max(0, RS**2 - (yS-y2)**2))
  p20 = (xLo, y2)
  p21 = (xB0 - dxB2, y2)
  p22 = (xB0 + dxB2, y2)
  p23 = (xB1 - dxB2, y2)
  p24 = (xB1 + dxB2, y2)
  p25 = (xS + dxS2, y2)

  y3 = yLo + 3*wd
  dxB3 = sqrt(max(0, RB**2 - (yB-y3)**2))
  dxS3 = sqrt(max(0, RS**2 - (yS-y3)**2))
  p30 = (xLo, y3)
  p31 = (xB0 - dxB3, y3)
  p32 = (xB0 + dxB3, y3)
  p33 = (xB1 - dxB3, y3)
  p34 = (xB1 + dxB3, y3)
  p35 = (xS + dxS3, y3)

  y4 = yLo + 4*wd
  dxS4 = 0
  p40 = (xLo, y4)
  p41 = (xS + dxS4, y4)
  
  # Joints of links near contours:

  ya = yLo + 0.5*wd
  dxSa = sqrt(max(0, RS**2 - (yS-ya)**2))
  pa0 = (xS + dxSa, ya)
  
  yb = yLo + 1.5*wd
  dxSb = sqrt(max(0, RS**2 - (yS-yb)**2))
  pb0 = (xS + dxSb, yb)
  
  yc = yLo + 2.5*wd
  dxBc = sqrt(max(0, RB**2 - (yB-yc)**2))
  dxSc = sqrt(max(0, RS**2 - (yS-yc)**2))
  pc0 = (xB0 - dxBc, yc)
  pc1 = (xB0 + dxBc, yc)
  pc2 = (xB1 - dxBc, yc)
  pc3 = (xB1 + dxBc, yc)
  pc4 = (xS + dxSc, yc)
  
  yd = yLo + 3.5*wd
  dxSd = sqrt(max(0, RS**2 - (yS-yd)**2))
  pd0 = (xS + dxSd, yd)

  tag = "f"
  
  # Raster traces:
  
  mv00 = move.make(p00, p01, mp_fill); move.set_name(mv00, "T00")
                                                         
  mv10 = move.make(p10, p11, mp_fill); move.set_name(mv10, "T10")
                                                         
  mv20 = move.make(p20, p21, mp_fill); move.set_name(mv20, "T20")
  mv21 = move.make(p22, p23, mp_fill); move.set_name(mv21, "T21")
  mv22 = move.make(p24, p25, mp_fill); move.set_name(mv22, "T23")
                                                         
  mv30 = move.make(p30, p31, mp_fill); move.set_name(mv30, "T30")
  mv31 = move.make(p32, p33, mp_fill); move.set_name(mv31, "T31")
  mv32 = move.make(p34, p35, mp_fill); move.set_name(mv32, "T33")
                                                         
  mv40 = move.make(p40, p41, mp_fill); move.set_name(mv40, "T40")

  # Link traces:

  mk00 = move.make(p00, p10, mp_fill); move.set_name(mk00, "L00")
  
  mk01 = move.make(p01, pa0, mp_fill); move.set_name(mk01, "L01")
  mka1 = move.make(pa0, p11, mp_fill); move.set_name(mka1, "La1")

  mk10 = move.make(p10, p22, mp_fill); move.set_name(mk10, "L10")
  
  mk11 = move.make(p11, pb0, mp_fill); move.set_name(mk11, "L11")
  mkb1 = move.make(pb0, p25, mp_fill); move.set_name(mkb1, "Lb1")

  mk20 = move.make(p20, p30, mp_fill); move.set_name(mk20, "L20")

  mk21 = move.make(p21, pc0, mp_fill); move.set_name(mk21, "L21")
  mkc1 = move.make(pc0, p31, mp_fill); move.set_name(mkc1, "Lc1")

  mk22 = move.make(p22, pc1, mp_fill); move.set_name(mk22, "L22")
  mkc2 = move.make(pc1, p32, mp_fill); move.set_name(mkc2, "Lc2")

  mk23 = move.make(p23, pc2, mp_fill); move.set_name(mk23, "L23")
  mkc3 = move.make(pc2, p33, mp_fill); move.set_name(mkc3, "Lc3")

  mk24 = move.make(p24, pc3, mp_fill); move.set_name(mk24, "L24")
  mkc4 = move.make(pc3, p34, mp_fill); move.set_name(mkc4, "Lc4")

  mk25 = move.make(p25, pc4, mp_fill); move.set_name(mk25, "L25")
  mkc5 = move.make(pc4, p35, mp_fill); move.set_name(mkc5, "Lc5")

  mk30 = move.make(p30, p40, mp_fill); move.set_name(mk30, "L30")

  mk31 = move.make(p35, pd0, mp_fill); move.set_name(mk31, "L31")
  mkd1 = move.make(pd0, p41, mp_fill); move.set_name(mkd1, "Ld1")

  # Join them into paths:
  PFS = []
  PLS = []
  
  if version == 0:
    # Each filling path is a single move:
    for mv in mv00, mv10, mv20,mv21,mv22, mv30,mv31,mv32, mv40:
      ph = path.from_moves([mv])
      path.set_name(ph, move.get_name(mv).replace("M","F"))
      PFS.append(ph)
    # Trivial path at starting point:
    phorg = path.from_points((org,), None, None)
    PFS.append(phorg)
    # Links may have 1 or 2 traces:
    for MLS in (
        (mk00,), (mk01,mka1,), (mk10,), (mk11,mkb1,),
        (mk20,), (mk21,mkc1,), (mk22,mkc2,), (mk23,mkc3,), (mk24,mkc4,), (mk25,mkc5,),
        (mk30,), (mk31,mkd1,) 
      ):
      ph = path.from_moves(MLS)
      path.set_name(ph, "L%d" % len(PLS))
      PLS.append(ph)
  else:
    # Single filling path -- define the list {TRS} of traces to be used, in order:
    if version == 1:
      # Scanline order with consistent directions, no links:
      TRS = [ mv00, mv10, mv20,mv21,mv22, mv30,mv31,mv32, mv40, ]
    elif version == 2:
      # Scanline order with alternating directions, links:
      TRS = [ 
        mv00, mk01,mka1, move.rev(mv10), 
        mv20,mv21,mv22, mk25,mkc5, 
        move.rev(mv32),move.rev(mv31),move.rev(mv30), mk30, mv40,
      ]
    elif version == 3:
      # Optimal(?) path:
      TRS = [
        mv00, mk01,mka1, move.rev(mv10), 
        move.rev(mv21), mk22,mkc2, mv31,
        mv22, mk25,mkc5, move.rev(mv32),
        move.rev(mv40), move.rev(mk30),
        mv30,
        move.rev(mkc1), move.rev(mk21),
        move.rev(mv20),
      ]
    
    # Now join them with jumps as needed:
    MVS = []; p = org;
    for mv in TRS:
      q = move.pini(mv)
      if p != q:
        jm = move.make(p, q, mp_jump)
        move.set_name(jm, move.get_name(mv).replace("M","J"))
        MVS.append(jm)
      MVS.append(mv)
      p = move.pfin(mv)
    ph = path.from_moves(MVS)
    path.set_name(ph, "F")
    PFS.append(ph)
  return PFS, PLS
  # ----------------------------------------------------------------------

def raster_figure(version, mp_cont, mp_fill, mp_jump):
  # Key dimensions and coordinates:
  wdc = move_parms.width(mp_cont)
  wdf = move_parms.width(mp_fill)

  off = (wdf + wdc)/2   # Offset of contour axis from filling axes/endpoints.

  # All coords refer to the axes and endpoints traces.

  xLo = 0                 # {X} coord of W side of filling.
  yLo = 0                 # {Y} coord of S side of filling.

  RS = 2*wdf              # Radius of  E semicircle of filling.
  xS = xLo + 9*wdf        # {X} coord of center of E semicirle
  yS = yLo + RS           # {Y} coord of center of E semicircle.

  yHi = yS + RS           # {Y} coord of N side.
 
  yB = yS + 0.5*wdf       # {Y} coord of centers of holes.
  xB0 = xLo + 3*wdf       # {X} coord of center of W hole.
  xB1 = xLo + 7*wdf       # {X} coord of center of E hole.
  RB = 1.5*wdf            # Radius of holes at the filling.

  xD = xLo + 4*wdf        # {X} coord of E edge of dent at SE corner.
  yD = yS                 # {Y} coord of N edge of dent at SE corner.

  PCS = raster_figure_contour(mp_cont, xLo,yLo, xS,yS,RS, yHi, xB0,xB1,yB,RB, xD,yD, off)
  PFS,PLS = raster_figure_filling(mp_fill, mp_jump, xLo,yLo, xS,yS,RS, yHi, xB0,xB1,yB,RB, xD,yD, wdf, version)
  return PCS,PFS,PLS
  # ----------------------------------------------------------------------

# MISCELLANEOUS RASTER FILLS FOR TESTING

def rasters_A(mp_fill, xdir, ydir, yphase, eps):

  # Endpoints in the standard coord system, unit step, zero phase, no perturbation:

  p00 = (0, 0); q00 = (3, 0)
  p20 = (5, 0); q20 = (9, 0)
  
  p01 = (0, 1); q01 = (9, 1)
  
  p02 = (1, 2); q02 = (2, 2)
  p22 = (4, 2); q22 = (6, 2)
  p42 = (8, 2); q42 = (9, 2)

  p03 = (1, 3); q03 = (2, 3)
  p23 = (4, 3); q23 = (6, 3)
  
  p = [ 
    [ p00, p20, p01, p02, p22, p42, p03, p23, ],
    [ q00, q20, q01, q02, q22, q42, q03, q23, ],
  ]
  
  # Scale, add phase, and rotate the points to the {xdir,ydir} system:
  wd = move_parms.width(mp_fill)
  OPHS = []
  for k in range(len(p[0])):
    prk = [None, None] # Scaled, rotated, perturbed endpoints
    for iend in range(2):
      pke = rn.scale(wd, p[iend][k])
      ke = 2*k + iend
      dpke = ( eps*sin(17*ke + 3), yphase + eps*cos(31*ke +5) ) # Perturbation to apply.
      ppke = rn.add(pke, dpke)
      prk[iend] = rn.mix(ppke[0], xdir, ppke[1], ydir)
    mvk = move.make(prk[0], prk[1], mp_fill);
    phk = path.from_moves((mvk,))
    OPHS.append(phk)

  return OPHS
  # ----------------------------------------------------------------------