function phi = basis_mn_tet4 ( t, n, p )

%*****************************************************************************80
%
%% basis_mn_tet4(): all bases at N points for a TET4 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the vertices of a tetrahedron.
%
%    It works directly with these coordinates, and does not refer to a
%    reference element.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Olgierd Zienkiewicz,
%    The Finite Element Method,
%    Sixth Edition,
%    Butterworth-Heinemann, 2005,
%    ISBN: 0750663200,
%    LC: TA640.2.Z54.
%
%  Input:
%
%    real T(3,4), the coordinates of the vertices.
%
%    integer N, the number of evaluation points.
%
%    real P(3,N), the points where the basis functions
%    are to be evaluated.
%
%  Output:
%
%    real PHI(4,N), the value of the basis functions
%    at the evaluation points.
%
  phi(1:4,1:n) = 0.0;
%
%           | x1 x2 x3 x4 |
%  Volume = | y1 y2 y3 y4 |
%           | z1 z2 z3 z4 |
%           |  1  1  1  1 |
%
  volume =                             ...
      t(1,1) * (                       ...
        t(2,2) * ( t(3,3) - t(3,4) )   ...
      - t(2,3) * ( t(3,2) - t(3,4) )   ...
      + t(2,4) * ( t(3,2) - t(3,3) ) ) ...
    - t(1,2) * (                       ...
        t(2,1) * ( t(3,3) - t(3,4) )   ...
      - t(2,3) * ( t(3,1) - t(3,4) )   ...
      + t(2,4) * ( t(3,1) - t(3,3) ) ) ...
    + t(1,3) * (                       ...
        t(2,1) * ( t(3,2) - t(3,4) )   ...
      - t(2,2) * ( t(3,1) - t(3,4) )   ...
      + t(2,4) * ( t(3,1) - t(3,2) ) ) ...
    - t(1,4) * (                       ...
        t(2,1) * ( t(3,2) - t(3,3) )   ...
      - t(2,2) * ( t(3,1) - t(3,3) )   ...
      + t(2,3) * ( t(3,1) - t(3,2) ) );

  if ( volume == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BASIS_MN_TET3 - Fatal error!\n' );
    fprintf ( 1, '  Element has zero volume.\n' );
    error ( 'BASIS_MN_TET4 - Fatal error!\n' );
  end
%
%             | xp x2 x3 x4 |
%  Phi(1,P) = | yp y2 y3 y4 | / volume
%             | zp z2 z3 z4 |
%             |  1  1  1  1 |
%
  phi(1,1:n) = (                           ...
      p(1,1:n) * (                         ...
        t(2,2)   * ( t(3,3)   - t(3,4) )   ...
      - t(2,3)   * ( t(3,2)   - t(3,4) )   ...
      + t(2,4)   * ( t(3,2)   - t(3,3) ) ) ...
    - t(1,2) * (                           ...
        p(2,1:n) * ( t(3,3)   - t(3,4) )   ...
      - t(2,3)   * ( p(3,1:n) - t(3,4) )   ...
      + t(2,4)   * ( p(3,1:n) - t(3,3) ) ) ...
    + t(1,3) * (                           ...
        p(2,1:n) * ( t(3,2)   - t(3,4) )   ...
      - t(2,2)   * ( p(3,1:n) - t(3,4) )   ...
      + t(2,4)   * ( p(3,1:n) - t(3,2) ) ) ...
    - t(1,4) * (                           ...
        p(2,1:n) * ( t(3,2)   - t(3,3) )   ...
      - t(2,2)   * ( p(3,1:n) - t(3,3) )   ...
      + t(2,3)   * ( p(3,1:n) - t(3,2) ) ) ) / volume;
%
%             | x1 xp x3 x4 |
%  Phi(2,P) = | y1 yp y3 y4 | / volume
%             | z1 zp z3 z4 |
%             |  1  1  1  1 |
%
  phi(2,1:n) = (                             ...
      t(1,1) * (                             ...
        p(2,1:n) * ( t(3,3)   - t(3,4) )     ...
      - t(2,3)   * ( p(3,1:n) - t(3,4) )     ...
      + t(2,4)   * ( p(3,1:n) - t(3,3) ) )   ...
    - p(1,1:n)   * (                         ...
        t(2,1)   * ( t(3,3)   - t(3,4) )     ...
      - t(2,3)   * ( t(3,1)   - t(3,4) )     ...
      + t(2,4)   * ( t(3,1)   - t(3,3) ) )   ...
    + t(1,3) * (                             ...
        t(2,1)   * ( p(3,1:n) - t(3,4) )     ...
      - p(2,1:n) * ( t(3,1)   - t(3,4) )     ...
      + t(2,4)   * ( t(3,1)   - p(3,1:n) ) ) ...
    - t(1,4) * (                             ...
        t(2,1)   * ( p(3,1:n) - t(3,3) )     ...
      - p(2,1:n) * ( t(3,1)   - t(3,3) )     ...
      + t(2,3)   * ( t(3,1)   - p(3,1:n) ) ) ) / volume;
%
%             | x1 x2 xp x4 |
%  Phi(3,P) = | y1 y2 yp y4 | / volume
%             | z1 z2 zp z4 |
%             |  1  1  1  1 |
%
  phi(3,1:n) = (                             ...
      t(1,1) * (                             ...
        t(2,2)   * ( p(3,1:n) - t(3,4) )     ...
      - p(2,1:n) * ( t(3,2)   - t(3,4) )     ...
      + t(2,4)   * ( t(3,2)   - p(3,1:n) ) ) ...
    - t(1,2) * (                             ...
        t(2,1)   * ( p(3,1:n) - t(3,4) )     ...
      - p(2,1:n) * ( t(3,1)   - t(3,4) )     ...
      + t(2,4)   * ( t(3,1)   - p(3,1:n) ) ) ...
    + p(1,1:n) * (                           ...
        t(2,1)   * ( t(3,2)   - t(3,4) )     ...
      - t(2,2)   * ( t(3,1)   - t(3,4) )     ...
      + t(2,4)   * ( t(3,1)   - t(3,2) ) )   ...
    - t(1,4) * (                             ...
        t(2,1)   * ( t(3,2)   - p(3,1:n) )   ...
      - t(2,2)   * ( t(3,1)   - p(3,1:n) )   ...
      + p(2,1:n) * ( t(3,1)   - t(3,2) ) ) ) / volume;
%
%             | x1 x2 x3 xp |
%  Phi(4,P) = | y1 y2 y3 yp | / volume
%             | z1 z2 z3 zp |
%             |  1  1  1  1 |
%
  phi(4,1:n) = (                             ...
      t(1,1) * (                             ...
        t(2,2)   * ( t(3,3)   - p(3,1:n) )   ...
      - t(2,3)   * ( t(3,2)   - p(3,1:n) )   ...
      + p(2,1:n) * ( t(3,2)   - t(3,3) ) )   ...
    - t(1,2) * (                             ...
        t(2,1)   * ( t(3,3)   - p(3,1:n) )   ...
      - t(2,3)   * ( t(3,1)   - p(3,1:n) )   ...
      + p(2,1:n) * ( t(3,1)   - t(3,3) ) )   ...
    + t(1,3) * (                             ...
        t(2,1)   * ( t(3,2)   - p(3,1:n) )   ...
      - t(2,2)   * ( t(3,1)   - p(3,1:n) )   ...
      + p(2,1:n) * ( t(3,1)   - t(3,2) ) )   ...
    - p(1,1:n) * (                           ...
        t(2,1)   * ( t(3,2)   - t(3,3) )     ...
      - t(2,2)   * ( t(3,1)   - t(3,3) )     ...
      + t(2,3)   * ( t(3,1)   - t(3,2) ) ) ) / volume;

  return
end
