function q = simpqual ( p, t, type )

%*****************************************************************************80
%
%% SIMPQUAL computes the simplex quality.
%
%  Discussion:
%
%    If TYPE = 1, the quality measure used computes the ratio of
%    the radii of the inscribed and circumscribed circles or spheres.
%
%    Only dimensions 1, 2 and 3 can be handled.  In 1D, the quality is always 1.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson. 
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,ND), the coordinates of a set of nodes.
%
%    Input, integer T(NT,1:ND+1), a list of the nodes which make up each 
%    simplex in the mesh.
%
%    Input, integer TYPE, specifies the quality measure:
%    1: Compute the radius ratio (default)
%    2: Use an approximate formula.
%
%    Output, real Q(NT), the simplex quality measure.
%
  if ( nargin < 3 )
    type = 1;
  end

  switch type
%
%  RADIUS RATIO
%
    case 1

      switch size ( p, 2 )

        case 1

          q = ones ( 1, size(t,2) );

        case 2

          a=sqrt(sum((p(t(:,2),:)-p(t(:,1),:)).^2,2));
          b=sqrt(sum((p(t(:,3),:)-p(t(:,1),:)).^2,2));
          c=sqrt(sum((p(t(:,3),:)-p(t(:,2),:)).^2,2));
          r=1/2*sqrt((b+c-a).*(c+a-b).*(a+b-c)./(a+b+c));
          R=a.*b.*c./sqrt((a+b+c).*(b+c-a).*(c+a-b).*(a+b-c));
          q = 2 * r ./ R;

        case 3

          d12=p(t(:,2),:)-p(t(:,1),:);
          d13=p(t(:,3),:)-p(t(:,1),:);
          d14=p(t(:,4),:)-p(t(:,1),:);
          d23=p(t(:,3),:)-p(t(:,2),:);
          d24=p(t(:,4),:)-p(t(:,2),:);
          d34=p(t(:,4),:)-p(t(:,3),:);
          v=abs(dot(cross(d12,d13,2),d14,2))/6;
          s1=sqrt(sum(cross(d12,d13,2).^2,2))/2;
          s2=sqrt(sum(cross(d12,d14,2).^2,2))/2;
          s3=sqrt(sum(cross(d13,d14,2).^2,2))/2;
          s4=sqrt(sum(cross(d23,d24,2).^2,2))/2;
          p1=sqrt(sum(d12.^2,2)).*sqrt(sum(d34.^2,2));
          p2=sqrt(sum(d23.^2,2)).*sqrt(sum(d14.^2,2));
          p3=sqrt(sum(d13.^2,2)).*sqrt(sum(d24.^2,2));
          q=216*v.^2./(s1+s2+s3+s4)./sqrt((p1+p2+p3).*(p1+p2-p3).* ...
                                    (p1+p3-p2).*(p2+p3-p1));

        otherwise

          error ( 'SIMPQUAL - Dimension not implemented.' );
    end
%
%  APPROXIMATE FORMULA
%
    case 2

     switch size(p,2)

       case 1

         q = ones ( 1, size(t,2) );

       case 2

         d12 = sum((p(t(:,2),:)-p(t(:,1),:)).^2,2);
         d13 = sum((p(t(:,3),:)-p(t(:,1),:)).^2,2);
         d23 = sum((p(t(:,3),:)-p(t(:,2),:)).^2,2);
         q = 4*sqrt(3)*abs(simpvol(p,t))./(d12+d13+d23);

       case 3

         d12 = sum((p(t(:,2),:)-p(t(:,1),:)).^2,2);
         d13 = sum((p(t(:,3),:)-p(t(:,1),:)).^2,2);
         d14 = sum((p(t(:,4),:)-p(t(:,1),:)).^2,2);
         d23 = sum((p(t(:,3),:)-p(t(:,2),:)).^2,2);
         d24 = sum((p(t(:,4),:)-p(t(:,2),:)).^2,2);
         d34 = sum((p(t(:,4),:)-p(t(:,3),:)).^2,2);
         q = 216*abs(simpvol(p,t))/sqrt(3)./(d12+d13+d14+d23+d24+d34).^(3/2);

       otherwise

         error ( 'SIMPQUAL - Dimension not implemented.' );
    end

  otherwise

    error ( 'SIMPQUAL - Incorrect type.' );

  end

  return
end