// Last edited on 2013-02-02 21:01:21 by stolfilocal

// Defines macros that compute a parametrizable vector field
// of a particular form (a translation term added to a periodic
// array of vortices) and the various terms of its Taylor expansion
// about the origin.

// ----------------------------------------------------------------------
// AFFINE (LINEAR + CONSTANT) FIELD 

#macro field_LIN_max(C,J, term)
  // Estimates the max magnitude of the field assuming the argument has 
  // length 1.
  #local ax = array[3];
  #local ax[0] = x;
  #local ax[1] = y;
  #local ax[2] = z;
  #if ((term = 0) | (term = 1))
    #local vmax_c = vlength(< C[0], C[1], C[2] >);
  #else
    #local vmax_c = < 0, 0, 0 >;
  #end
  
  // Find the row of the linear part with maximum norm:
  #local vmax_r = 0;
  #local i = 0;
  #while (i < 3)
    // Get row {i} of the relevant matrix:
    #switch (term)
    #case (0)
    #case (1)
    #case (4)
    #case (9)
      // Linear term is {J} itself:
      #local row = < J[i][0], J[i][1], J[i][2] >;
      #break
    #case (2)
    #case (3)
      // Linear term is zero:
      #local row = < 0, 0, 0 >;
      #break
    #case (5)
      // Linear term is the antisymmetric part of {J}, namely {(J-J')/2}:
      #local row = 0.5 * < J[i][0]-J[0][i], J[i][1]-J[1][i], J[i][2]-J[2][i] >;
      #break
    #case (6)
      // Linear term is the symmetric part of {J}, namely {(J+J')/2}:
      #local row = 0.5 * < J[i][0]+J[0][i], J[i][1]+J[1][i], J[i][2]+J[2][i] >;
      #break
    #case (7)
      // Linear term is the symmetric isotropic component of {J}, namely {trace(J)/3*I}:
      #local row = (J[0][0]+J[1][1]+J[2][2])/3 * ax[i];
      #break
    #case (8)
      // Linear term is the traceles symmetric component of {J} namely {(J+J')/2 - trace(J)/3*I}:
      #local row = 0.5 * < J[i][0]+J[0][i], J[i][1]+J[1][i], J[i][2]+J[2][i] > - (J[0][0]+J[1][1]+J[2][2])/3 * ax[i];
      #break
    #else
      #error "invalid Taylor {term} code"
    #end
    #local vrow = vlength(row); 
    
    #if (vrow > vmax_r) 
      #local vmax_r = vrow;
    #end
    #local i = i + 1;
  #end
  vmax_c + vmax_r;
#end

#macro field_LIN(C,J, xp,yp,zp)
  // The official steady-state flow velocity field.
  // Returns the velocity vector at point {xp,yp} 
  // and time {tt} in a static flow that has the following
  // first-order components at the origin
  //   * mean velocity {VX,VY}
  //   * rotation with angular velocity {W}
  //   * expansion with relative rate {D}
  //   * shear with diagonal {S,-S} and off-diagonal 0
  //   * extra quadratic term {Q(x^2+y^2,x^2-y^2)}
  // The velocity gradient is mostly constant near the origin 
  // but changes far away.  The parameter RM determines the
  // distance when this happens.

  #local vp = field_LIN_term(C,J, 0, xp,yp,zp)
  vp
#end

// ----------------------------------------------------------------------
// TAYLOR DECOMPOSITION
// The following terms define the Taylor terms of the flow 
// about the origin:

#macro field_LIN_term(C,J, term, xp,yp,zp)
  #local p = array[3];
  #local p[0] = xp;
  #local p[1] = yp;
  #local p[2] = zp;
  #local c = array[3];
  #local c[0] = C[0];
  #local c[1] = C[1];
  #local c[2] = C[2];
  #switch (term)
  #case (0)
    // There are no non-linear terms:
    #local vp = field_LIN_taylor_affine(c,J, p);
    #break
  #case (1)
    #local vp = field_LIN_taylor_affine(c,J, p);
    #break
  #case (2)
    // There are no non-linear terms:
    #local vp = < 0, 0, 0 >;
    #break
  #case (3)
    #local vp = field_LIN_taylor_constant(c,J, p);
    #break
  #case (4)
    #local vp = field_LIN_taylor_linear(c,J, p);
    #break
  #case (5)
    #local vp = field_LIN_taylor_rotation(c,J, p);
    #break
  #case (6)
    #local vp = field_LIN_taylor_strain(c,J, p);
    #break
  #case (7)
    #local vp = field_LIN_taylor_expansion(c,J, p);
    #break
  #case (8)
    #local vp = field_LIN_taylor_shear(c,J, p);
    #break
  #case (9)
    // There are no non-linear terms:
    #local vp = field_LIN_taylor_linear(c,J, p);
    #break
  #else
    #error "invalid Taylor {term} code"
  #end
  vp
#end

#macro field_LIN_taylor_affine(c,J, p)
  // The Taylor approximation of the flow field at the origin.
  // The constant term plus the linear term
  #local vv = array[3];
  #local i = 0;
  #while (i < 3)
    #local vv[i] = J[i][0]*p[0] + J[i][1]*p[1] + J[i][2]*p[2] + c[i];
    #local i = i + 1;
  #end
  < vv[0], vv[1], vv[2] >
#end

#macro field_LIN_taylor_constant(c,J, p)
  // The constant term at the origin.
  
  #local vv = < c[0], c[1], c[2] >;
  < vv[0], vv[1], vv[2] >
#end

#macro field_LIN_taylor_linear(c,J, p)
  // The linear term at the origin
  #local vv = array[3];
  #local i = 0;
  #while (i < 3)
    #local vv[i] = J[i][0]*p[0] + J[i][1]*p[1] + J[i][2]*p[2];
    #local i = i + 1;
  #end
  < vv[0], vv[1], vv[2] >
#end

#macro field_LIN_taylor_rotation(c,J, p)
  // The rotation (antisymmetric) part of the linear term at the origin
  
  // Rotation term:
  #local vv = array[3];
  #local i = 0;
  #while (i < 3)
    #local vv[i] = 0.5*((J[i][0] - J[0][i])*p[0] + (J[i][1] - J[1][i])*p[1] + (J[i][2] - J[2][i])*p[2]);
    #local i = i + 1;
  #end
  < vv[0], vv[1], vv[2] >
#end

#macro field_LIN_taylor_strain(c,J, p)
  // The strain (symmetric) part of the linear term at the origin
  #local vv = array[3];
  #local i = 0;
  #while (i < 3)
    #local vv[i] = 0.5*((J[i][0] + J[0][i])*p[0] + (J[i][1] + J[1][i])*p[1] + (J[i][2] + J[2][i])*p[2]);
    #local i = i + 1;
  #end
  < vv[0], vv[1], vv[2] >
#end

#macro field_LIN_taylor_expansion(c,J, p)
  // The expansion (trace) part of the linear term at the origin
  #local sm = 0;
  #local i = 0;
  #while (i < 3)
    #local sm = sm + J[i][i];
    #local i = i + 1;
  #end
  #local D = sm/3.0;
  < D*p[0], D*p[1], D*p[2] >
#end

#macro field_LIN_taylor_shear(c,J, p)
  // The shear (symmetric traceless) part of the linear term at the origin
  #local vp = 
    field_LIN_taylor_strain(c,J, p) -
    field_LIN_taylor_expansion(c,J, p);
  vp
#end