      FUNCTION CAngle (CNum1)
      INCLUDE 'STDHDR.FOR'
      COMPLEX CNum1
      REAL NearZero
      PARAMETER (NearZero = 1E-20)


      IF (ABS(REAL(CNum1)) .LT. NearZero) THEN
        IF (ABS(IMAG(CNum1)) .LT. NearZero)  CAngle = 0
        IF (IMAG(CNum1) .LT. -(NearZero))  CAngle = 3.0 * pi / 2.0
        IF (IMAG(CNum1) .GT. NearZero) CAngle = pi / 2.0
      ELSE
        IF (ABS(IMAG(CNum1)) .LT. NearZero) THEN
          IF (REAL(CNum1) .LT. -(NearZero))  CAngle = pi
          IF (REAL(CNum1) .GT. NearZero)  CAngle = 0
        ELSE
          IF (REAL(CNum1) .GT. NearZero .AND.
     +        IMAG(CNum1) .GT. NearZero)
     +           CAngle = ATAN(IMAG(CNum1) / REAL(CNum1))
          IF (REAL(CNum1) .LT. -NearZero .AND.
     +        IMAG(CNum1) .GT. NearZero)
     +           CAngle = pi - ATAN(-IMAG(CNum1) / REAL(CNum1))
          IF (REAL(CNum1) .LT. -NearZero .AND.
     +        IMAG(CNum1) .LT. -NearZero)
     +           CAngle = pi + ATAN(IMAG(CNum1) / REAL(CNum1))
          IF (REAL(CNum1) .GT. NearZero .AND.
     +        IMAG(CNum1) .LT. -NearZero)
     +         CAngle = 2.0 * pi - ATAN(-IMAG(CNum1) / REAL(CNum1))
        END IF
      END IF
      END !FUNCTION CAngle

      SUBROUTINE CSwap (n1, n2)
      REAL n1, n2 ,temp
        temp = n1
        n1 =  n2
        n2 = temp
      END !SUB CSwap


      FUNCTION CMag (CNum1)
      COMPLEX CNum1
        CMag = SQRT((REAL(CNum1) ** 2) + (IMAG(CNum1) ** 2))
      END !FUNCTION CMag



      SUBROUTINE CMatAdd (m1, m2, numrow, numcol, m3)
      INCLUDE 'STDHDR.FOR'
      COMPLEX m1(0:maxr,0:maxc), m2(0:maxr, 0:maxc), m3(0:maxr, 0:maxc)
      INTEGER numrow,numcol
      INTEGER i, j

      DO i = 0, numrow - 1
        DO j = 0, numcol - 1
          CALL ComplexMath(m1(i, j), 0, m2(i, j), m3(i, j))
        END DO
      END DO
      END !SUB CMatAdd



      SUBROUTINE ComplexMatInvert (matdata, numcol, det, invary)
      INCLUDE 'STDHDR.FOR'
      REAL matdata(0:maxr,0:maxc),  invary(0:maxr, 0:maxc)
      INTEGER numcol
      REAL det

      INTEGER i, j, k, l, lerow, lecol, l1, iErr
      REAL piv, t, leval
      INTEGER pivlst[ALLOCATABLE](:,:)
      LOGICAL pivchk[ALLOCATABLE](:)
      ALLOCATE(pivlst(0:maxc,0:1),STAT=iErr)
      CALL CheckMem(iErr)
      ALLOCATE(pivchk(0:maxc),STAT=iErr)
      CALL CheckMem(iErr)

      det = 1.0

      DO i = 0, numcol - 1
        pivchk(i) = .false.
        DO j = 0, numcol - 1
          invary(i, j) = matdata(i, j)
        END DO
      END DO

      DO i = 0, numcol - 1
        leval = 0.0
        DO j = 0, numcol - 1
           IF (.NOT. pivchk(j)) THEN
             DO k = 0, numcol - 1
               IF (.NOT. pivchk(k)) THEN
                 IF (ABS(invary(j, k)) .GT. leval) THEN
                   lerow = j
                   lecol = k
                   leval = ABS(invary(j, k))
                 END IF
               END IF
             END DO
           END IF
         END DO

         pivchk(lecol) = .TRUE.
         pivlst(i, 0) = lerow
         pivlst(i, 1) = lecol

        IF (lerow .NE. lecol) THEN
           det = -det
           DO l = 0, numcol - 1
             CALL CSwap (invary(lerow, l), invary(lecol, l))
           END DO
         END IF

         piv = invary(lecol, lecol)
         det = det * piv
         IF (det .GT. 1E30) THEN
           det = 1E30       ! {prevents numeric overflow}
         END IF
         invary(lecol, lecol) = 1.0
         DO l = 0, numcol - 1
           invary(lecol, l) = invary(lecol, l) / piv
         END DO
         DO l1 = 0, numcol - 1
           IF (l1 .NE. lecol)  THEN
             t = invary(l1, lecol)
             invary(l1, lecol) = 0
             DO l = 0, numcol - 1
               invary(l1, l) = invary(l1, l) - invary(lecol, l) * t
             END DO
           END IF
         END DO
       END DO

       DO i = 0, numcol - 1
         l = numcol - i - 1
         IF (pivlst(l, 0) .NE. pivlst(l, 1)) THEN
           lerow = pivlst(l, 0)
           lecol = pivlst(l, 1)
           DO k = 0, numcol - 1
             CALL CSwap( invary(k, lerow), invary(k, lecol))
           END DO
         END IF
       END DO
       DEALLOCATE(pivlst,STAT=iErr)
       CALL CheckDealloc(iErr)
       DEALLOCATE(pivchk,STAT=iErr)
       CALL CheckDealloc(iErr)
      END !SUB





      SUBROUTINE CMatInvert (matdata, numcol, det, invary)
      INCLUDE 'STDHDR.FOR'
      COMPLEX matdata(0:maxr,0:maxc), invary(0:maxr, 0:maxc)
      REAL det
      REAL coef[ALLOCATABLE](:,:), inv[ALLOCATABLE](:,:)
      INTEGER numcol, i,j,nn,ii,jj,iErr
      ALLOCATE (coef(0:maxr, 0:maxc),inv(0:maxr, 0:maxc), STAT=iErr)

      DO i = 0, numcol - 1
        ii = 2*i
        DO j = 0, numcol - 1
          jj = 2*j
          coef(ii, jj) = REAL(matdata(i, j))
          coef(ii, jj+ 1) = -IMAG(matdata(i, j))
          coef(ii + 1, jj ) = IMAG(matdata(i, j))
          coef(ii + 1, jj + 1) = REAL(matdata(i, j))
         END DO
      END DO
      nn = 2 * numcol

      CALL ComplexMatInvert(coef, nn, det, inv)

      DO i = 0, numcol - 1
        ii = 2*i
        DO j = 0, numcol - 1
          jj = 2*j
          invary(i, j) = CMPLX(inv(ii, jj),inv(ii + 1, jj))
        END DO
      END DO
      DEALLOCATE (coef,inv, STAT=iErr)
      END !SUB CMatInvert


      SUBROUTINE CMatPrint (m1, m, n, matname)
      INCLUDE 'STDHDR.FOR'
      COMPLEX m1(0:maxr,0:maxc)
      INTEGER m,n
      CHARACTER * 80 matname
      INTEGER i, j

      WRITE (*,*) matname
      WRITE (*,*)
      DO i = 0, m - 1
        DO j = 0, n - 1
          WRITE (*,432) '[',REAL(m1(i,j)), IMAG(m1(i,j)),'j]'
        END DO
        WRITE (*,*)
      END DO
      WRITE (*,*)
432    FORMAT (1X, A1, F7.2, 2X, F7.2, A2 \)
      END !SUBROUTINE



      SUBROUTINE CMatProd (m1, m2, l, m, n, m3)
      INCLUDE 'STDHDR.FOR'
      COMPLEX m1(0:maxr,0:maxc), m2(0:maxr, 0:maxc), m3(0:maxr, 0:maxc)
      INTEGER l,m,n
      COMPLEX temp,temp2
      INTEGER i,j,k

      DO i = 0, l - 1
        DO j = 0, n - 1
          m3(i,j) = CMPLX(0,0)
          DO k = 0, m - 1
            CALL ComplexMath(m1(i, k), 2, m2(k, j), temp)
            CALL ComplexMath(m3(i, j), 0, temp, temp2)
            m3(i,j) = temp2
          END DO
        END DO
      END DO
      END !SUBROUTINE CMatProd



      SUBROUTINE CMatScalarProd (m1, sval, numrow, numcol, m2)
      INCLUDE 'STDHDR.FOR'
      COMPLEX m1(0:maxr,0:maxc), sval, m2(0:maxr, 0:maxc)
      INTEGER i, j

      DO i = 0, numrow - 1
        DO j = 0, numcol - 1
          CALL ComplexMath(sval, 2, m1(i, j), m2(i, j))
        END DO
      END DO
      END !SUB CMatScalarProd



      SUBROUTINE CMatTranspose (m1, numrow, numcol, m2)
      INCLUDE 'STDHDR.FOR'
      COMPLEX m1(0:maxr,0:maxc),  m2(0:maxr, 0:maxc)
      INTEGER   i, j

      DO i = 0, numrow - 1
        DO j = 0, numcol - 1
          m2(j, i) = m1(i, j)
        END DO
      END DO
      END !SUBROUTINE CMatTranspose





      SUBROUTINE ComplexMath (c1, op, c2, result)
      COMPLEX c1,c2,result
      INTEGER op
      REAL tempr, tempi

      SELECT CASE (op)
         CASE (0)
           result = c1 + c2
         CASE (1)
           result = c1 - c2
         CASE (2)
           tempr = REAL(c1) * REAL(c2) - IMAG(c1) * IMAG(c2)
           tempi = REAL(c1) * IMAG(c2) + IMAG(c1) * REAL(c2)
           result = CMPLX( tempr, tempi)
         CASE (3)
           tempr = (REAL(c1)*REAL(c2)+IMAG(c1)*IMAG(c2))/
     +                    ((REAL(c2) ** 2) + (IMAG(c2) ** 2))
           tempi = (IMAG(c1) * REAL(c2) - REAL(c1) * IMAG(c2))/
     +                    ((REAL(c2) ** 2) + (IMAG(c2) ** 2))
           result = CMPLX( tempr, tempi)
        END SELECT

      END !SUBROUTINE ComplexMath


