function [R, Mz] = mfs_warpR(x, y, z, evec)

# usage: [R, Mz] = mfs_warpR(x, y, z, evec)
#
# Input  x(4), y(4), z(4)  Nodal point coordinates in element
#                          coordinate system
#        evec(3, 4)        Vectors along element edges
#
# Output R(8)              Matrix to compute force R
#        Mz(8)             Matrix to compute modification of Mz
#
# The function computes matrix R to compute the forces and moments
# needed to make moment vectors parallel to the tangential plane.
#
# --------------------------------------------------------------------

  i1 = [4, 1, 2, 3];  % Indices of left edge vectors

  v  = cross(evec(:, i1), evec);
  vp = v(1 : 2, :) ./ v(3, :);

  Mz(2 : 2 : 8) = vp(2, :); Mz(1 : 2 : 7) = vp(1, :);

  s = x(2) + x(3) - x(1) - x(4) + y(3) + y(4) - y(1) - y(2);
  R = Mz / s;

  Mz = -Mz;

end
