#! /usr/bin/gawk -f
# Last edited on 1999-02-13 20:33:10 by stolfi

BEGIN {
  abort = -1;
  usage = ( \
    "ring-place.awk \\\n" \
    "  -v iRadius=NUM \\\n" \
    "  -v oRadius=NUM \\\n" \
    "  -v spacing=NUM \\\n" \
    "  > ring.coords" \
  );
  
  # Computes pebble positions and colors. Writes for each pebble in
  # the ring a line in the format
  #
  #   X Y RSZ GUN 
  #
  # where RSZ is the relative pebble size in [0 _ 1] and
  # GUN is the color gun index (0..2).
  
  pi = 3.1415926;
  twopi = 2*pi;
  
  if (iRadius !~ /^[0-9]*[.]?[0-9]+$/)      { arg_error("bad or missing \"iRadius\""); }
  if (oRadius !~ /^[0-9]*[.]?[0-9]+$/)      { arg_error("bad or missing \"oRadius\""); }
  if (spacing !~ /^[0-9]*[.]?[0-9]+$/)      { arg_error("bad or missing \"spacing\""); }
  
  compute_particles();
  
}

function compute_particles(  \
  nr,npb,dt,dta,i,r,rr,rho,rsz,spc,per,nt,nta,j,tau,jj,t,tt,x,y, \
  c0,c1,c2,x0,x1,x2,y0,y1,y2,gun \
)
{
  printf "computing particle positions...\n" > "/dev/stderr";
  split("", ocolor); split("", oposx); split("", oposy);
  nr = int((oRadius - iRadius)/(spacing * sqrt(3)/2));
  npb = 0;
  nt = 0;
  for (i = 0; i < nr; i ++)
    { r = iRadius + (oRadius - iRadius)*(i+0.5)/nr;
      split("", color);
      
      # relative size:
      rho = sqrt((i+0.5)/nr);
      rsz = 0.10 + 0.10*rho + 0.80*(4*rho*(1-rho));
      
      # adjusted spacing:
      spc = (1.40 - (1-rho)*(1-rho))*spacing;

      per = 2 * 3.1415926 * r;
      nta = nt; nt = 3*int(per/spc/3);
      dt = dt + 0.25*(rand()-0.5) + int(nt*rand());
      while (dt < 0) { dt += nt; }
      while (dt >= nt) { dt -= nt; }
      printf "i = %d nt = %d dt = %6.4f\n", i, nt, dt > "/dev/stderr"
      gun = 0;
      for (j = 0; j < nt; j ++)
        { 
          # Nominal azimuth (in turns).
          # Try to offset it by 1/2 step rel. to previous row.
          t = (j + 0.5*i + dt + 0.25*(rand()-0.5))/nt;
          
          # Compute "gun" (color) and adjust coordinates x,y:
          if (i == 0) 
            { gun = j % 3;
              rr = r + (rand()-0.5)*spc;
              x = rr * cos(tt*twopi);
              y = rr * sin(tt*twopi);
            }
          else
            { # Look up coordinates and colors of 
              # nearest neighbors in previous layer:
              jj = int(nta*(t - 0.50001/nt));
              while (jj < 0) { jj += nta; } jj %= nta;
              c0 = ocolor[jj]; x0 = oposx[jj]; y0 = oposy[jj];
              
              jj = int(nta*(t + 0.50001/nt));
              while (jj < 0) { jj += nta; } jj %= nta;
              c1 = ocolor[jj]; x1 = oposx[jj]; y1 = oposy[jj];
              
              # Pick the nominal position along the radius:
              x = (x0 + x1)/2; y = (y0 + y1)/2;
              rr = sqrt(x*x + y*y) + spc + 0.5*(rand()-0.5)*spc;
              x = rr * cos(t*twopi);
              y = rr * sin(t*twopi);
              
              # Adjust x and y:
              adjust_coordinates(x,y,x0,y0,x1,y1,spc);
              x = adj_x; y = adj_y;
              
              # Compute a distinctive color for the current particle:
              if (c0 == c1) 
                { if (c0 == gun) 
                    { c2 = (gun + 1) % 3; }
                  else
                    { c2 = 3 - c0 - gun; }
                }
              else
                { c2 = 3 - c0 - c1; }
              if ((c2 < 0) || (c2 >= 3))
                { prog_error(i,j,("color computation failed: c0=" c0 " c1=" c1 " c2=" c2)); }
              gun = c2;
            }
            
          # Save color and coordinates for next iteration.
          jj = int(nt*t);
          while (jj < 0) { jj += nt; } jj %= nt;
          color[jj] = gun; posx[jj] = x; posy[jj] = y;
          # make sure there are no holes:
          jj++; jj %= nt;
          if ((jj > 0) || (! (jj in color))) 
            { color[jj] = gun; posx[jj] = x; posy[jj] = y; }

          printf "%6.3f %6.3f  %5.3f %d\n", x, y, rsz, gun
          npb++;
        }
      # Prepare for next iteration:
      
      split("", ocolor); split("", oposx); split("", oposy);
      for (jj=0; jj < nt; jj++)
        { ocolor[jj] = color[jj]; oposx[jj] = posx[jj]; oposy[jj] = posy[jj]; }
    }
  printf "%d pebbles generated\n", npb > "/dev/stderr";
}

function adjust_coordinates(x,y,x0,y0,x1,y1,spc,    dx,dy,xc,yc,d,d2,h)
{
  # Returns in "adj_x, adj_y" the nominal pebble coordinates "(x,y)",
  # adjusted by taking into account the neighbors "(x0, y0)" and
  # (x1,y1)" in the previous ring. They will not deviate by more than
  # half the step "spc".
  dx = x1-x0; dy = y1-y0;
  d2 = dx*dx + dy*dy;
  if (d2 < 0.000001*spc)
    { # push the point (x,y) away from (x0,y0)
      dx = x-x0; dy = y-y0;
      if ((dx == 0) && (dy == 0)) { dx = 0.00001; }
      d = sqrt(dx*dx + dy*dy);
      xc = x0 + spc*dx/d;
      yc = y0 + spc*dy/d;
    }
  else
    { # make an isosceles triangle with area = spc^2*sqrt(3)/4
      h = (spc*spc/d2)*(sqrt(3)/2);
      xc = x0 + dx/2 + dy*h;
      yc = y0 + dy/2 - dx*h;
    }

  # Add a small random perturbation (less than spc/4):
  dx = 1; dy = 1;
  while (dx*dx + dy*dy > 0.25) { dx = rand()-0.5; dy = rand()-0.5; }
  xc += 0.5*dx*spc;
  yc += 0.5*dy*spc;
  
  # Now prevent it from wandering too far off the current position:
  dx = xc - x;
  dy = yc - y;
  d = sqrt(dx*dx + dy*dy);
  h = (spc/2)/(d + (spc/2));
  adj_x = x + h*dx;
  adj_y = y + h*dy;
}


function prog_error(i,j,msg)
{
  printf "*** i=%d j=%d: %s\n", i, j, msg > "/dev/stderr";
  abort = 1; exit abort;
}

function arg_error(msg)
{
  printf "*** %s\n", msg > "/dev/stderr";
  abort = 1; exit abort;
}