function [ d, ifault ] = detq ( a, n )

%*****************************************************************************80
%
%% detq() computes the determinant of an orthogonal matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 January 2020
%
%  Author:
%
%    Original FORTRAN77 version by J C Gower.
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    J C Gower,
%    Algorithm AS 82:
%    The determinant of an orthogonal matrix,
%    Applied Statistics,
%    Volume 24, Number 1, 1975, page 150-153.
%
%  Input:
%
%    real A(N,N), the orthogonal matrix stored by rows or columns.
%
%    integer N, the order of the matrix.
%
%  Output:
%
%    real D, the determinant of A.
%
%    integer IFAULT, 
%    0, no error occurred.
%    1, an error was detected.
%
  ifault = 0;
  d = 0.0;

  tol = 0.0001;

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'detq(): Fatal error!\n' );
    fprintf ( 1, '  n <= 0\n' );
    ifault = 1;
    return
  end

  a2 = zeros ( n * n, 1 );
  k = 0;
  for i = 1 : n
    for j = 1 : n
      k = k + 1;
      a2(k) = a(i,j);
    end
  end

  d = 1.0;
  r = 1;

  for k = 2 : n + 1

    q = r;
    x = a2(r);
    y = sign ( x );
    d = d * y;
    y = - 1.0 / ( x + y );
    x = abs ( x ) - 1.0;

    if ( tol < abs ( x ) )

      if ( 0.0 < x )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'detq(): Fatal error!\n' );
        fprintf ( 1, '  x < 0.0\n' );
        fprintf ( 1, '  x = %g\n', x );
        ifault = 1;
        return
      end

      if ( k == n + 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'detq - Fatal error!\n' );
        fprintf ( 1, '  k == n + 1\n' );
        ifault = 1;
        return
      end

      for i = k : n
        q = q + n;
        x = a2(q) * y;
        p = r;
        s = q;
        for j = k : n
          p = p + 1;
          s = s + 1;
          a2(s) = a2(s) + x * a2(p);
        end

      end

    end

    r = r + n + 1;

  end

  return
end

