function Q = mfs_warpQ(x, y, z)

# usage: Q = mfs_warpQ(x, y, z)
#
# Input  x(4), y(4), z(4)  Nodal point coordinates in element
#                          coordinate system
# Output Q(2, 4)           Matrix to compute Q forces
#
# The function computes matrix Q that relates the vertical warping
# correction forces to the nodal point forces.
#
# --------------------------------------------------------------------

  C = [y(1) - y(3), y(2) - y(4); x(1) - x(3), x(2) - x(4)];
  Z(1, 2 : 2 : 8) = z'; 
  Z(2, 1 : 2 : 7) = z';
  Q = C \ Z;

end
