#! /usr/bin/gawk -f
# Last edited on 1999-02-13 22:07:07 by stolfi

BEGIN {
  abort = -1;
  usage = ( \
    "ring-twiddle.awk \\\n" \
    "    -v pokerX=NUM \\\n" \
    "    -v pokerY=NUM \\\n" \
    "  < ring-locs.txt \\\n" \
    "  > ring-props.txt" \
  );
  
  # Given the X,Y positions of each pebble, computes the 
  # height Z, the slopes SX = dZ/dX and SY = dZ/dY, and
  # the "age" of the pebble in the wave..
  #
  # For each pebble, reads a line in the format X Y RSZ GUN,
  # and writes a line X Y Z SX SY RSZ AGE GUN
  
  pi = 3.1415926; twopi = 2*pi;

  # Debugging parameters
  image = 0;
  plot = 0;

  # Wave parameters:
  
  bowSize = 5.0;
  bowCurv = -0.20;
  bowAngl = 0.80;
  tipSize = 0.07;
  waveLen = 1.0;
  waveAmp = 0.8;
  # rippling = 0.075;
  rippling = 0.25

  # these parameters are dimensionless
  waveTan = 0.05

  if (pokerX   !~ /^[-+]?[0-9]*[.]?[0-9]+$/) { arg_error("bad or missing \"pokerX\""); }
  if (pokerY   !~ /^[-+]?[0-9]*[.]?[0-9]+$/) { arg_error("bad or missing \"pokerY\""); }
  
  if (plot) { make_test_plot(); exit 0; }
  if (image) { make_test_image(); exit 0; }
  printf "computing particle height and tilt...\n" > "/dev/stderr";
}

function make_test_plot(   nx,ny,sx,sy,step,ix,iy,x,y,z,age)
{
  # Write a file suitable for gnuplot
  
  nx = 100;
  ny = 100;
  sx = 3*bowSize / nx;
  sy = 3*bowSize / ny;
  step = (sx > sy ? sx : sy);
  printf "writing test plot...\n" > "/dev/stderr";
  
  for (ix = 0; ix < nx; ix ++)
    { for (iy = 0; iy < ny; iy ++)
        { 
          # Cartesian coordinates: x, y
          x = pokerX + (ix - 0.1*nx)*step;
          y = pokerY + (iy - 0.5*ny)*step;

          compute_wave_args(x,y);
          z = compute_wave_z();
          age = compute_wave_age();
          
          printf "%7.5f %7.5f %7.5f  %7.5f %7.5f  %7.5f %7.5f %7.5f  %7.5f\n", 
            x, y, z,  wave_u, wave_v,  wave_g, wave_h, wave_m, age > "debug.plot"
        }
      printf "\n" > "debug.plot"
    }
}  

function make_test_image(   image,nx,ny,sx,sy,step,x,y,ix,iy,z,age)
{
  # Create and write a debugging image:
  
  split("", image);
  nx = 256;
  ny = 256;
  sx = 2*bowSize / nx;
  sy = 2*bowSize / ny;
  step = (sx > sy ? sx : sy);
  printf "computing test image...\n" > "/dev/stderr";
  
  for (ix = 0; ix < nx; ix ++)
    { for (iy = 0; iy < ny; iy ++)
        { 
          # Cartesian coordinates: x, y
          x = pokerX + (ix - 0.1*nx)*step;
          y = pokerY + (iy - 0.5*ny)*step;

          compute_wave_args(x,y);
          z = compute_wave_z()
          age = compute_wave_age()
          
          image[iy,ix] = age
        }
    }
  dump_image(image, nx,ny)
}  

function dump_image(im, nx, ny,   ix,iy,v,maxv)
{
  printf "normalizing image...\n" > "/dev/stderr";
  maxv = 0.000001;
  for (iy = ny-1; iy >= 0; iy--)
    { for (ix = 0; ix < nx; ix++)
        { v = im[iy,ix] + 0;
          v = (v > 0 ? v : -v);
          if (v > maxv) { maxv = v; }
        }
    }
  printf "maxv = %f\n", maxv > "/dev/stderr";
  printf "writing image to debug.pgm...\n" > "/dev/stderr";
  printf "P2\n" > "debug.pgm";
  printf "%d %d\n", nx, ny > "debug.pgm";
  printf "255\n" > "debug.pgm";
  for (iy = ny-1; iy >= 0; iy--)
    { for (ix = 0; ix < nx; ix++)
        { v = int(127.5*(im[iy,ix]/maxv + 1));
          if (v < 0) { v = 0; }
          if (v > 255) { v = 255; }
          printf " %d", v > "debug.pgm";
        }
      printf "\n" > "debug.pgm"; 
    }
  close("debug.pgm");
  printf "done.\n" > "/dev/stderr";
}

(abort >= 0) { exit abort;} 

/./ { 
  if (NF != 4) { file_error("bad format"); }
  x = $1;
  y = $2;
  rsz = $3;
  gun = $4;

  compute_wave_args(x,y);
  z = compute_wave_z();
  age = compute_wave_age();
  
  # if (z != 0) 
  #   { printf "z = %6.4f wave_t = %6.4f age = %6.4f\n", \
  #       z, wave_t, age > "/dev/stderr";
  #   }

  dx = 0.0001;
  compute_wave_args(x+dx, y);
  sx = (compute_wave_z() - z)/dx;

  dy = 0.0001;
  compute_wave_args(x, y+dy);
  sy = (compute_wave_z() - z)/dy;

  printf "%6.3f %6.3f %6.4f  %5.3f %5.3f  %5.3f %5.3f %d\n", 
    x,y,z, sx, sy, rsz, age, gun;
  next;
}

function compute_wave_args(x,y,  xx,yy,ca,sa,psi,q2max,umin,vmax,scale,p2,q2,s2,f)
{
  # Computes the wave paramaters for point (x,y):
  #   "wave_x_raw" = re-centered and rotated coordinate, absolute
  #   "wave_y_raw" = re-centered and rotated coordinate, absolute
  #   "wave_x" = re-centered, rotated and rectified longitudinal coordinate, absolute
  #   "wave_y" = re-centered, rotated and rectified latitudinal coordinate, absolute
  #   "wave_u" = scaled and translated wave_x (u=1 critical pt)
  #   "wave_v" = scaled wave_y
  #   "wave_h" = relative shock amplitude (0__1)
  #   "wave_g" = relative ripple amplitude (0__1)  
  #   "wave_a" = distance from upper reference line (rel tipSize)
  #   "wave_b" = distance from lower reference line (rel tipSize)
  #   "wave_t" = time since shock front passed through point
  #   "wave_m" = position indicator (0 = front/outside, 1=axis).
  #   "wave_r" = smoothed distance from focus (rel tipSize, 1__+oo )
  # from  the cartesian coordinates (x,y)
  
  # Uses the conformal wave coordinates: p (along), q (across)
  # where q = const is a family of hyperbolas with foci (±1,0)
  # and p = const is a family of ellipses with same foci.
  
  # Subtract hand position:
  x = x - pokerX;
  y = y - pokerY;
  
  # Rotate coordinate system:
  ca = cos(bowAngl); sa = sin(bowAngl);
  wave_x_raw = x * ca + y * sa;
  wave_y_raw = y * ca - x * sa;
  x = wave_x_raw; y = wave_y_raw;
  
  # Remove curvature of wave axis
  psi = pi*x/bowSize
  yy = y - bowCurv*bowSize*sin(psi);
  xx = x + pi*cos(psi)*atan2(bowCurv*yy,bowSize);
  wave_x = xx; wave_y = yy;

  if ((xx == 0) && (yy == 0)) { xx = 0.0001; }

  # These parameters are all relative to tipSize:
  wave_a = (waveTan*xx - yy)/tipSize;
  wave_b = (waveTan*xx + yy)/tipSize;
  wave_r = sqrt(xx*xx + yy*yy + tipSize*tipSize)/tipSize;

  # Define the "wave scale" coordinates (u,v) = (wave_u,wave_v).
  # These are a uniform scaling and horizontal shifting of (xx,yy).
  # so that the "focal point" of the wave u=1 is at xx = 0,
  # and the tip u=umin is at xx = -tipSize.

  # The outer envelope of the wave is defined 
  # by the equation wave_q2(u,v) = waveTan^2,
  # which is approximately |v/u| = waveTan for large fixed u.
  # The leftmost point is at v = 0 and u = sqrt(1/(1 + waveTan^2)).
  # Therefore
  
  q2max = waveTan*waveTan;
  umin = sqrt(1/(1 + q2max));
  scale = tipSize/(1 - umin);

  wave_u = xx/scale + 1;
  wave_v = yy/scale;
  
  # The local "wave clock" wave_t is the time since the wavefront
  # passed the current point. It is computed assuming that the
  # wavefront moves towards negative u, at such a speed that the
  # distance from tip to critical point is traversed in 1 unit of
  # time. Negative times mean that the wave hasn't reached the 
  # current point yet.
  #
  # wave_ufront(v,q2max) is the u coordinate of the point on the
  # wavefront with the current v, i.e. such that wave_q2(u,v) = q2max
  
  wave_t = (wave_u - wave_ufront(wave_v, q2max))/(1-umin);
  
  # The parametres g, h, and m are functions of the 
  # ratio s = q/qmax between the q coordinate and front's q = qmax.
  # They are such that g^2 + h^2 = m^2/f(r) where m is the 
  # indicator function sqrt(1 - s^4) and f is a decay function
  # of the radius r.
  
  wave_h = 0; 
  wave_g = 0;
  wave_m = 0;
  
  # Crude clip to wave boundary:
  vmax = waveTan*wave_u;
  if ((wave_u <= umin) || (wave_v > vmax) || (wave_v < -vmax)) { return; }
  # p2 = wave_p2(wave_u,wave_v);
  q2 = wave_q2(wave_u,wave_v);

  if (q2 >= q2max) { return; }
  
  s2 = q2/q2max 
  if (s2 > 1) 
    { wave_g = 0; wave_h = 0; wave_m = 0; }
  else 
    { wave_g = sqrt((1 + s2 - s2*s2)*(1 - s2)); 
      wave_h = s2*sqrt(1 - s2);
      wave_m = sqrt(1 - s2*s2);
    }
  
  # Introduce the 1/f(r) decay. Note that dec is 1 at the critical
  # point, and is about 2/k at points k*tipSize away from it,
  # for large k.
  
  f = sqrt(15 + wave_r*wave_r)/4
  wave_h /= f;
  wave_g /= f;
}

function compute_wave_z_test(  x,y,r,res)
{
  x = twopi*wave_x_raw/waveLen;
  y = twopi*wave_y_raw/waveLen;
  r = sqrt(x*x + y*y);
  if (r == 0) 
    { res = 1; }
  else
    { res = sin(r)/r; }
  return(waveAmp * res);
}

function compute_wave_z(   s,ea,eb,ec,wa,wb,wc,aa,ab,ac)
{
  # wave_z is the relative height of the wave (-A__+A where A < 2)
  #
  # Must call compute_wave_args first

  s = twopi/waveLen;
  
  ea = s*(0.800*wave_x + 0.600*wave_y);
  eb = s*(0.921*wave_x - 0.400*wave_y);
  ec = s*(0.644*wave_x - 0.765*wave_y);
  
  aa = +0.80*rippling;
  ab = -0.50*rippling;
  ac = -0.33*rippling;

  wa = aa*cos(0.667*ea);
  wb = ab*cos(1.000*eb);
  wc = ac*cos(1.500*ec);

  return waveAmp * (wave_h*sqrt(1-aa*aa-ab*ab-ac*ac) + wave_g*(wa + wb + wc));
}

function wave_p2(u,v,  u2,v2,A,B,C,D,S,P)
{
  # Distance along wave
  u2 = u*u;
  v2 = v*v;
  A = 1;
  B = u2 + v2 - 1;
  C = -u2;
  P = qsolve(A,B,C);
  return P
}

function wave_q2(u,v,  u2,v2,A,B,C,D,S,Q)
{
  # This parameter is zero on the ray { (u,0) : u >= 1},
  # grows approximately like (v/u)^2 for large fixed u,
  # grows to infinity as u tends to zero,
  # and is (1 - u^2)/u^2 for v = 0 and u between 0 and 1.
  
  u2 = u*u;
  v2 = v*v;
  if (v2 == 0) 
    { if (u2 > 1) { Q = 0; } else { Q = (1-u2)/u2 } }
  else
    { A = u2; B = u2 - v2 - 1; C = -v2;
      Q = qsolve(A,B,C);
    }
  return Q;
}

function wave_ufront(v,q2,  v2,q4)
{
  # Computes u such that wave_q2(u,v) = q2
  q4 = q2*q2;
  v2 = v*v;
  return sqrt((v2 + q2 + v2*q2)/(q2 + q4));
}

function compute_wave_age(  t,ep,em)
{
  return wave_m*wave_m;
}

function tanh(t,  s,ep)
{ 
  # Hyperbolic tangent of t.
  # result is in [-1 __ +1] as t ranges in [-oo __ +oo].
  
  if (t < 0) { t = -t; s = -1; } else { s = +1; }
  if (t > 15)
    { return s; }
  else 
    { ep = exp(-2*t);
      return s*(1 - ep)/(1 + ep);
    }
}

function qsolve(A,B,C, BB,AA,AC,M,S)
{
  # Largest root of A*x^2 + B*x + C = 0
  BB = B*B;  AA = A*A;  AC = A*C;
  M = (AC < 0 ? -AC : AC); 
  M = (M > AA ? M : AA);
  if (M < 0.00000001*BB)
    { return(-C/B); }
  else if (B == 0)
    { if (A > 0) { return (sqrt(-C/A)); } else { return (-sqrt(-C/A)); } }
  else if (C == 0)
    { if (B < 0) { return (-B/A); } else { return(0); } }
  else
    { S = (B - sqrt(BB - 4*AC))/2; return (-S/A); }
}

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

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