c/* %W%      latest revision %G% %U% */
      subroutine cmatin(a, b, n, ndim)
*******************************************************************************

c
c
c *** finds the inverse of complex matrix h=(rea,aima)=(a,b)
c
c *** code provided by norman bardsley to kwon+tabakin s prog bopit
c *** n=no of elements in a assigned(and in the sub=dimension)
c *** matrix a,b(real+imag parts) destroyed by call
c *** the other arguments are dummy just to save space

      implicit real*8 (a-h, o-z)

      dimension  indexa(200), indexb(200), ipivot(200)
      dimension  a(ndim,ndim), b(ndim,ndim)

      equivalence (irow,jrow), (icolum,jcolum)

c ***
c *** detr=1.0
c *** rhl add to define tempi for real potential

      tempi = 0.

c *** deti=0.0

      do 10 j = 1,n
  10     ipivot(j) = 0

      do 70 i = 1,n
         xmax = 0.0
         do 35 j = 1,n
            if (ipivot(j) - 1) 15, 35, 15

  15        do 30 k=1,n
               if (ipivot(k) - 1) 20, 30, 90
  20           if (abs(a(j,k)) .lt. 1.0e-20) go to 30
               xjk = abs(a(j,k)) + abs(b(j,k))
               if (xmax - xjk)  25,30,30
  25           irow = j
               icolum = k
               xmax = xjk
  30        continue

  35  continue

      ipivot(icolum) = ipivot(icolum) + 1
      if (irow .eq. icolum) go to 50

c *** detr=-deti
c *** deti=-deti

      do 45 l = 1,n
         swap = a(irow,l)
         swapi = b(irow,l)
         a(irow,l) = a(icolum,l)
         b(irow,l) = b(icolum,l)
         b(icolum,l) = swapi
  45     a(icolum,l) = swap

  50  indexa(i) = irow
      indexb(i) = icolum
      pivotr = a(icolum,icolum)
      pivoti = b(icolum,icolum)

c *** temp=detr*pivotr-deti*pivoti
c *** deti=detr*pivoti+deti*pivotr
c *** detr=temp

      a(icolum,icolum) = 1.0
      b(icolum,icolum) = 0.0
      if (pivoti .eq. 0.0) tempr = 1.0/pivotr
      if (pivoti .ne. 0.0) then
         tempr =  pivotr/(pivotr * pivotr + pivoti * pivoti)
         tempi = -pivoti/(pivotr * pivotr + pivoti * pivoti)
      endif

      do 55 l = 1,n
         temp = a(icolum,l) * tempr - b(icolum,l) * tempi
         b(icolum,l) = a(icolum,l) * tempi + b(icolum,l) * tempr
  55     a(icolum,l) = temp

      do 70 l1 = 1,n
         if (l1 - icolum) 60, 70, 60
  60     tempa = a(l1,icolum)
         tempb = b(l1,icolum)
         a(l1,icolum) = 0.0
         b(l1,icolum) = 0.0

         do 65 l = 1,n
            b(l1,l) = b(l1,l) - a(icolum,l) * tempb
     $                - b(icolum,l) * tempa
            a(l1,l) = a(l1,l) - a(icolum,l) * tempa
     $                + b(icolum,l) * tempb
  65     continue
  70  continue

      do 85 i = 1,n
         l = n + 1 - i
         if (indexa(l) - indexb(l)) 75, 85, 75
  75     jrow = indexa(l)
         jcolum = indexb(l)

         do 80 k = 1,n
            swap = a(k,jrow)
            swapi = b(k,jrow)
            a(k,jrow) = a(k,jcolum)
            b(k,jrow) = b(k,jcolum)
            a(k,jcolum) = swap
            b(k,jcolum) = swapi
  80     continue
  85  continue

  90  return
      end