function ref = triangle_order6_physical_to_reference ( t, n, phy )

%*****************************************************************************80
%
%% triangle_order6_physical_to_reference() maps a physical point to a reference point.
%
%  Discussion:
%
%    Given the vertices of an order 6 physical triangle and a point
%    (X,Y) in the physical triangle, the routine computes the value
%    of the corresponding image point (R,S) in reference space.
%
%    The mapping from (R,S) to (X,Y) has the form:
%
%      X(R,S) = A1 * R * R + B1 * R * S + C1 * S * S
%             + D1 * R     + E1 * S     + F1
%
%      Y(R,S) = A2 * R * R + B2 * R * S + C2 * S * S
%             + D2 * R     + E2 * S     + F2
%
%  Reference Element T3:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  6  5
%    |  |   \
%    |  |    \
%    0  1--4--2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real T(2,6), the coordinates of the vertices.  The vertices are 
%    assumed to be the images of (0,0), (1,0), (0,1), (1/2,0), (1/2,1/2)
%    and (0,1/2), in that order.
%
%    integer N, the number of points to transform.
%
%    real PHY(2,N), the coordinates of points in the physical space.
%
%  Output:
%
%    real REF(2,N), the coordinates of the corresponding
%    points in the reference space.
%

%
%  Set iteration parameters.
%
  it_tol = 0.000001;
  it_max = 10;

  for i = 1 : 2

    a(i) =   2.0 * t(i,1) + 2.0 * t(i,2)                ...
           - 4.0 * t(i,4);

    b(i) =   4.0 * t(i,1)                               ...
           - 4.0 * t(i,4) + 4.0 * t(i,5) - 4.0 * t(i,6);

    c(i) =   2.0 * t(i,1)                + 2.0 * t(i,3) ...
                                         - 4.0 * t(i,6);

    d(i) = - 3.0 * t(i,1) -       t(i,2)                ...
           + 4.0 * t(i,4);

    e(i) = - 3.0 * t(i,1)                -       t(i,3) ...
                                         + 4.0 * t(i,6);
    f(i) =         t(i,1);

  end
%
%  Initialize the points by inverting the linear map.
%
  ref(1:2,1:n) = triangle_order3_physical_to_reference ( t(1:2,1:3), n, phy(1:2,1:n) );
%
%  Carry out the Newton iteration.
%
  for j = 1 : n

    for it = 1 : it_max

      fun(1:2) = a(1:2) * ref(1,j) * ref(1,j) ...
               + b(1:2) * ref(1,j) * ref(2,j) ...
               + c(1:2) * ref(2,j) * ref(2,j) ...
               + d(1:2) * ref(1,j) ...
               + e(1:2) * ref(2,j) ...
               + f(1:2) ...
               - phy(1:2,j)';

      fun_norm = sqrt ( fun(1) * fun(1) + fun(2) * fun(2) );

      if ( fun_norm <= it_tol )
        break
      end

      jac(1:2,1) = 2 * a(1:2)' * ref(1,j) +     b(1:2)' * ref(2,j) + d(1:2)';
      jac(1:2,2) =     b(1:2)' * ref(1,j) + 2 * c(1:2)' * ref(2,j) + e(1:2)';


      dx(1:2) = - jac(1:2,1:2) \ fun(1:2)';

      ref(1:2,j) = ref(1:2,j) + dx(1:2)';

    end

  end

  return
end
