// Last edited on 2013-01-31 23:51:25 by stolfilocal

// The client must predefine the {field} macro.
// Defines macro {sampler_grid} that generates an image of 
// a term of the flow defined by the macro {field}.

#include "particle_with_arrow.inc"

#macro sampler_grid(zp,tt, term,  Rg,Rclip, vmag,Ra, dp,mj, frp,txp, fro, fsa,fra,txa)
  // Generates a scene consisting of a grid of particles with
  // arrows showing the direction and magnitude of the specified {term} of velocity field {field} at time {tt}.
  
  // {zp}    Z of plane to be sampled.
  // {tt}    Time to be sampled.
  // {term}  Selects the term of Taylor's expansion that will be shown.

  // {Rg}    Grid radius. The particle grid will span from {-Rg} to {+Rg} in {X} and {Y}.
  // {Rclip} If positive, clips the grid to a disk of radius {round}.
  
  // {vmag}  Max velocity magnitude to assume for scaling arrow lengths.
  // {Ra}    Max desirable arrow length for scaling arrow lengths.

  // {dp}    The particle spacing.
  // {mj}    Magnitude of jitter in particle positions.

  // {frp}   Extra factor for particle radius. Normally 1. 
  // {txp}   Texture for particles.

  // {fro}   Extra factor for particle and arrow radius at origin. Normally 1 or 0. 
  
  // {fsa}   Extra factor for arrows lengths. Normally 1. 
  // {fra}   Extra factor for arrows shaft radius. Normally 1. 
  // {txa}   Texture for velocity arrows.

  // There will be a particle at the origin iff {np} is odd and {fro > 0}.
  
  // The parameters {fsa,fra,frp} are extra scale factors
  // for the length and radius of arrow shafts and the radius of particles.
  // The parameter {fro} is an extra factor
  // for the radius of the ball and radius of the particle at origin (if any).
  // If these parameters are set to 1, a velocity vector with length {vmag} 
  // is mapped to an arrow whose length is half the particle spacing; and the 
  // arrow and ball radii should be appropriate.
  
  #local Rd = Rclip - 0.25*Ra; // Actual clipping radius, to avoid particles running into ring.
  
  #local dx = dp;            // Particle spacing in X direction.
  #local dy = dx*sqrt(3)/2;  // Spacing of particle rows.

  #local hnx = int(Rg/dx + 0.999999);  // Half particle count along X.
  #local hny = int(Rg/dy + 0.999999);  // Half row count along Y.
  
  // sampler_grid_test()
  
  union{

    #local iy = -hny;
    #while (iy <= +hny)
      #local mx = 0.5*mod(iy,2); // Relative row shift.
      #local ix = -hnx;
      #while (ix + mx <= +hnx)
        // Compute coord of point in {[-Rg _ +Rg]^2}: 
        #local xp = (ix+mx)*dx;
        #local yp = iy*dy;
        #local ev = sampler_grid_jitter(ix,iy);
        // #debug concat("  ", str(ev.x,9,6), " ", str(ev.y,9,6), " ", str(ev.z,9,6), " # ev \n") 
        #local pp = < xp, yp, zp > + < mj, mj, 0> * ev;

        #if ((Rclip = 0) | ((xp*xp + yp*yp) < Rd*Rd))
        
          #local vp = field(pp.x,pp.y,pp.z,tt, term);

          #if ((ix != 0) | (iy != 0))
            object{ particle_with_arrow(pp,vp, vmag,Ra, frp,txp, fsa,fra,txa) }
          #else
            #if (fro > 0)
              object{ particle_with_arrow(pp,vp, vmag,Ra, fro*frp,txp, fsa,fro*fra,txa) }
            #end
          #end
        #end

        #local ix = ix + 1;
      #end
      
      #local iy = iy + 1;
    #end

  }
#end

#macro sampler_grid_test()
  #local hn = 50;
  #local iy = -hn;
  #while (iy <= +hn)
    #local ix = -hn;
    #while (ix <= +hn)
      #local ev = sampler_grid_jitter(ix,iy);
      #debug concat("  ", str(ev.x,9,6), " ", str(ev.y,9,6), " ", str(ev.z,9,6), " # ev \n") 
      #local ix = ix + 1;
    #end
    #local iy = iy + 1;
  #end
#end

#macro sampler_grid_jitter(ix,iy)
  // A pseudorandom function of {ix,iy} to {\RR^3}. Each coord is independent Gaussian
  // with zero mean and approx unit deviation.
  
  #local ev = < 0, 0, 0 >;
  #local rr = 0;
  #local kk = 0;
  #local nt = 36;
  #while (kk < nt)
    #local rr = sin((100*ix*(2*ix-1) + 321*iy*(2*iy-1) + 1)*rr + 1);
    #local ev = < ev.y + rr, ev.y, ev.z >;
    #local rr = sin((210*ix*(2*ix-1) + 317*iy*(2*iy-1) + 1)*rr + 1);
    #local ev = < ev.x, ev.z + rr, ev.z >;
    #local rr = sin((151*ix*(2*ix-1) + 417*iy*(2*iy-1) + 1)*rr + 1);
    #local ev = < ev.x, ev.y, ev.x + rr >;
    #local kk = kk+1;
  #end
  (1.15 * ev / sqrt(nt))
#end