/*
diderotc --exec bug028.diderot && ./bug028
*/

input int sizeI = 90;
input int sizeJ = 80;
input real radius = 1.0;
input bool domain = false;
vec2 margin = [0.7,0.985];
input int which = 0;

function vec2 idx2wrl (vec2 idx)
  = [lerp(-radius, radius, 2, idx[0], sizeI-3),
     lerp(radius, -radius, 2, idx[1], sizeJ-3)];

/* bspln3 scaled to max out at 1.0 */
function real bspl (real x)
 = (0.0 if |x| > 2 else (8 - 12*|x| + 6*|x|^2 - |x|^3 if |x| > 1 else (4 - 6*|x|^2 + 3*|x|^3)))/4;

function real bump (real x) = (0.0 if x < -margin[1]
                               else (bspl(lerp(-2,0,-margin[1],x,-margin[0])) if x < -margin[0]
                               else (1.0 if x < margin[0]
                               else (bspl(lerp(0,2,margin[0],x,margin[1])) if x < margin[1]
                               else 0.0))));

function vec2 ff (vec2 pos) {
  real xx = pos[0];
  real yy = pos[1];
  vec2 V = [0,0];
  V = ([ xx, yy] if 0 == which else
      ([-xx,-yy] if 1 == which else
      ([-xx, yy] if 2 == which else
      ([ xx,-yy] if 3 == which else
      ([ yy, xx] if 4 == which else
      ([-yy,-xx] if 5 == which else
      ([-yy, xx] if 6 == which else
      ([ yy,-xx]))))))));
  return bump(xx)*bump(yy)*V;
}

vec2 oo = idx2wrl([0,0]);
vec2 uu = idx2wrl([1,0]) - idx2wrl([0,0]);
vec2 vv = idx2wrl([0,1]) - idx2wrl([0,0]);

// rows appear as contiguous vectors
tensor[2,3] i2w = [[uu[0], vv[0], oo[0]],
                   [uu[1], vv[1], oo[1]]];

strand sample (int ii, int jj) {
  vec3 posi = [ii, jj, 1];
// BUG: 2-by-3 times 3-vector multiply not implemented?
uncaught exception Fail [Fail: unsupported matrix-vector multiply]
  raised at common/phase-timer.sml:76.50-76.52
  raised at common/phase-timer.sml:76.50-76.52
  raised at c-util/tree-to-c.sml:185.30-185.71
*/
  vec2 posw = i2w•posi;
  output vec2 out = [0,0];

  update {
    if (ii == 0 && jj == 0) {
      //print("unu save -i out.nrrd -f nrrd -o out.nhdr\n");
      print("space dimension: 2\n");
      print("space directions: none (", uu[0], ",", uu[1], ") (", vv[0], ",", vv[1], ")\n");
      print("space origin: (", oo[0], ",", oo[1], ")\n");
    }
    out = posw if domain else ff(posw);
    stabilize;
  }
}

initially [sample(ii, jj) | jj in 0..sizeJ-1, ii in 0..sizeI-1];
//                           slow    --->  to  -->  fast
//                           in both strand creation order, and
//                           saving output raster order