#	function determ -- Calculates determinant of a square matrix
#
#   Version  1.1B   31-JAN-85 16:13:30	   BALL 	  Module designated a UTILITY
#		     1-SEP-92		R.L. WILLIAMSON    Modified for STSDAS use
#		     21-JAN 93		 R.L. Williamson   Converted to SPP
#   Description:
#   ------------
#   This REAL*8 function returns the value of the determinant of the
#   square matrix specified. It is based on the DETERM program to be
#   found in Bevington (appendix B, page 294). CAUTION! The original
#   array will be destroyed by this routine during diagonalization!!
#
#   Inputs:
#   -------
#   ARRAY
#   ACTUAL_ARRAY_SIZE
#   ORDER_OF_ARRAY
#
#   Outputs: (result of function)
#   --------
#   DETERMINANT_OF_ARRAY
#
#
#	Version    Date        Author	  Description
#	   6	26-NOV-1984  Z.G. Levay   Updated PDL
#	   5	 3-NOV-1983  PAT MURPHY   Submitted for inspection
#	   4	 2-NOV-1983  PAT MURPHY   Fixed bug
#	   3	 1-NOV-1983  PAT MURPHY   Submitted for inspection
#	   2	22-SEP-1983  PAT MURPHY   Add Fortran code
#	   1	22-SEP-1983  PAT MURPHY   Prolog Generation
#
#  .
#  DO for each matrix column
#	.
#	IF the element on the diagonal is Zero
#		.
#		IF the matrix is singular
#			.
#			Set the Determinant to Zero
#			Terminate
#			.
#		ENDIF	(Singular matrix?)
#		.
#		DO for each row in column
#			.
#			Swap array elements
#			.
#		ENDDO	(Rows)
#		.
#	ENDIF	(Zero diagonal element?)
#	.
#	Subtract the current row from lower rows to get a diagonal matrix
#	.
#  ENDDO
#  .
#  RETURN
#  END
#
#

double procedure determ (array, size, norder)

int norder, size, i, j, k, k1
bool singular
double array [size, size], save, aik, akk,dter
define termin_ 10

begin

#  Initial value
	dter = 1.

	do k = 1, norder {

#	  Swap columns if diagonal element is zero

	  if (array[k,k] == 0.) {
#  See if matrix is singular
	    singular = true
	    do j = k, norder {
	      singular = singular && array[k,j] == 0.
	    }
#  If any column is all zero
	    if (singular) {
#  then set result to zero
#  and quit now
            return (0.0d0)
	    }
	    do i = k, norder {
#  swap array elements
	      save	  = array [i,j]
	      array [i,j] = array [i,k]
	      array [i,k] = save
	    }
	    dter = - dter
	  }

#	Subtract k'th row from lower rows so we get diagonal matrix

	  dter = dter * array [k,k]
	  if (k < norder) {
	    akk = array [k,k]
	    k1 = k + 1
	    do i = k1, norder {
	      aik = array [i,k]
	      do j = k1, norder {
		array [i,j] = array [i,j] - aik * array [k,j] / akk
	      }
	    }
	  }
	}

return (dter)
end