// Computes the velocity vector at a given point 

#macro flow_velocity(crv,A,B,C,D, xp,yp)

  #if (crv = 0)
    // Straight pathlines:
    #local yo = yp - D*xp;
    #local vox = B + C*yo;
    #local voy = D*vox;
    #local vp = < vox, voy, 0 >;
  #else
    // Circular pathlines:
    #local R = 1/crv;
    #local dxp = xp;
    #local dyp = R-yp;
    #local Ro = sqrt(dxp*dxp + dyp*dyp);
    #local tp = degrees(atan2(dxp,dyp)); // [sic]
    #local yo = R - Ro;
    #local vox = A/Ro + B + C*yo;
    #local voy = 0;
    #local vp = vrotate(< vox, 0, 0 >, tp*z);
  #end
  vp

#end