// Last edited on 2013-01-20 10:50:19 by stolfilocal

// Defines macros that compute a parametrizable steady-state 2D flow
// 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.

// ----------------------------------------------------------------------
// FULL NON-LINEAR FIELD

#macro flow_vel_TSC_max(VX,VY,W,D,S,Q)
  // Estimates the max magnitude of the velocity.
  #local vmax_const = sqrt(VX*VX  + VY*VY);
  #local vmax_rot = W*0.5;
  #local vmax_exp = D+S;
  #local vmax = vmax_const + sqrt(vmax_rot*vmax_rot + vmax_exp*vmax_exp);
  (0.25*vmax)
#end

#macro flow_vel_TSC(VX,VY,W,D,S,Q, xp,yp)
  // 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.
  
  // Rotation term:
  #local vx_r = +W/pi*sin(pi*yp)*cos(pi*xp);
  #local vy_r = -W/pi*sin(pi*xp)*cos(pi*yp);
  
  // Strain term:
  #local vx_s = +(D+S)/pi*sin(pi*xp)*cos(pi*yp);
  #local vy_s = +(D-S)/pi*sin(pi*yp)*cos(pi*xp);
  
  // Extra quadratic term: 
  #local vx_q = Q*(xp*xp + yp*yp);
  #local vy_q = Q*(xp*xp - yp*yp);
  
  // Total velocity:
  #local vp = < VX + vx_r + vx_s + vx_q, VY + vy_r + vy_s + vy_q, 0 >;
  vp
#end

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

#macro flow_vel_TSC_term(VX,VY,W,D,S,Q, term, xp,yp)
  #switch (term)
  #case (0)
    #local vp = flow_vel_TSC(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (1)
    #local vp = flow_vel_TSC_taylor_affine(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (2)
    #local vp = flow_vel_TSC_taylor_residual(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (3)
    #local vp = flow_vel_TSC_taylor_constant(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (4)
    #local vp = flow_vel_TSC_taylor_linear(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (5)
    #local vp = flow_vel_TSC_taylor_rotation(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (6)
    #local vp = flow_vel_TSC_taylor_strain(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (7)
    #local vp = flow_vel_TSC_taylor_expansion(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (8)
    #local vp = flow_vel_TSC_taylor_shear(VX,VY,W,D,S,Q, xp,yp);
    #break
  #case (9)
    #local vp = flow_vel_TSC_taylor_relative(VX,VY,W,D,S,Q, xp,yp);
    #break
  #else
    #error "invalid Taylor {term} code"
  #end
  vp
#end

#macro flow_vel_TSC_taylor_relative(VX,VY,W,D,S,Q, xp,yp)
  // The Taylor approximation of the flow field at the origin.
  // The full flow minus the constant term
  
  #local vp = 
    flow_vel_TSC(VX,VY,W,D,S,Q, xp,yp) - 
    flow_vel_TSC_taylor_constant(VX,VY,W,D,S,Q, xp,yp);
  vp
#end

#macro flow_vel_TSC_taylor_affine(VX,VY,W,D,S,Q, xp,yp)
  // The Taylor approximation of the flow field at the origin.
  // The constant term plus the linear term
  
  #local vp = 
    flow_vel_TSC_taylor_constant(VX,VY,W,D,S,Q, xp,yp) + 
    flow_vel_TSC_taylor_linear(VX,VY,W,D,S,Q, xp,yp);
  vp
#end

#macro flow_vel_TSC_taylor_residual(VX,VY,W,D,S,Q, xp,yp)
  // The quadratic and higher terms of the Taylor
  // approximation of the flow field at the origin.
  
  #local vp = 
    flow_vel_TSC(VX,VY,W,D,S,Q, xp,yp) - 
    flow_vel_TSC_taylor_affine(VX,VY,W,D,S,Q, xp,yp);
  vp
#end

#macro flow_vel_TSC_taylor_constant(VX,VY,W,D,S,Q, xp,yp)
  // The constant term at the origin.
  
  #local vp = < VX, VY, 0 >;
  vp
#end

#macro flow_vel_TSC_taylor_linear(VX,VY,W,D,S,Q, xp,yp)
  // The linear term at the origin
  
  #local vp = 
    flow_vel_TSC_taylor_rotation(VX,VY,W,D,S,Q, xp,yp) + 
    flow_vel_TSC_taylor_strain(VX,VY,W,D,S,Q, xp,yp);
  vp
#end

#macro flow_vel_TSC_taylor_rotation(VX,VY,W,D,S,Q, xp,yp)
  // The rotation part of the linear term at the origin
  
  // Rotation term:
  #local vx_r = +W*yp;
  #local vy_r = -W*xp;
  
  // Total velocity:
  #local vp = < vx_r, vy_r, 0 >;
  vp
#end

#macro flow_vel_TSC_taylor_strain(VX,VY,W,D,S,Q, xp,yp)
  // The strain part of the linear term at the origin
  
  #local vp = 
    flow_vel_TSC_taylor_expansion(VX,VY,W,D,S,Q, xp,yp) + 
    flow_vel_TSC_taylor_shear(VX,VY,W,D,S,Q, xp,yp);
  vp
#end

#macro flow_vel_TSC_taylor_expansion(VX,VY,W,D,S,Q, xp,yp)
  // The expansion part of the linear term at the origin
  
  // Expansion term:
  #local vx_s = +D*xp;
  #local vy_s = +D*yp;
  
  // Total velocity:
  #local vp = < vx_s, vy_s, 0 >;
  vp
#end

#macro flow_vel_TSC_taylor_shear(VX,VY,W,D,S,Q, xp,yp)
  // The shear part of the linear term at the origin
  
  // Shear term:
  #local vx_s = +S*xp;
  #local vy_s = -S*yp;
  
  // Total velocity:
  #local vp = < vx_s, vy_s, 0 >;
  vp
#end