function [x, y, TE, T33] = mfs_esys(coor, opt)

# usage: [x, y, TE, T33] = mfs_esys(coor, opt)
#
# Input  coor(3, 3)   Coordinates of element nodes
#        opt          if 1, return 2 x 9 transformation matrix
#                     if 2, return 3 x 9 transformation matrix
#
# Output x(3)         x-coordinates of element nodes
#        y(3)         y-coordinates of element nodes
#        TE(:, 18)    Matrix to transform displacements to element
#                     coordinate system
#        T33(3, 3)    Elementary transformation matrix to element
#                     coordinate system
#
# The function computes the x- and y-coordinates of the element
# nodal points with respect to the element coordinate system and the
# transformation matrix from the global to the element coordinate
# system for triangular elements.
#
# --------------------------------------------------------------------

# Basis vectors of element coordinate system

  bx = coor(2, :) - coor(1, :);
  bx = bx / norm(bx);
  c2 = coor(3, :) - coor(1, :);
  bz = cross(bx, c2);
  bz = bz / norm(bz);
  by = cross(bz, bx);

  T33 = [bx; by; bz];

# Coordinates in element coordinate system

  relcoor = coor - coor(1, :);

  x = relcoor * bx';
  y = relcoor * by';

# Transformation matrix

  if (opt == 1)
     ix   = 1 : 2; incx = 2;
     TE1  = T33(1: 2, :);
     TE   = sparse(6, 9);
  else
     ix   = 1 : 3; incx = 3;
     TE1  = T33;
     TE    = sparse(9, 9);
  end 

  iy = 1 : 3;

  for n = 1 : 3
      TE(ix, iy) = TE1;
      ix += incx; iy += 3;
  end

end
