c/* @(#)matchc.f	1.1      latest revision 6/23/89 14:48:52 */
c// Lu, 7/16/91, changed this routine for direct matching of t matrix.
c// Lu, 12/8/92, rewrote this routine for coupled channel cases;.
c// Lu, 1/30/93, this routine is primarily written for L=J(+-)1, S=1 case,
c// however, it's also ture for L=J, S=0.5(+-)0.5 case.
c// Input: 3 complex t matrix; output: same 3 complex t matrix
c// BB convention
      subroutine matchc2(tr00,ti00,tr02,ti02,tr20,ti20,
     &	tr22,ti22,jdum,ldmax,aovera,xgam,ljtab)
      implicit real*8  (a-h,o-y)
      implicit complex*16 (z)
	  real*8 k0
	  real*8 Mnuc, Mnuc2, MN, Mp, Mp2
      real*8 fc(100),fcp(100),gc(100),gcp(100),scoul(100)
      real*8 fc0(100),fcp0(100),gc0(100),gcp0(100)
       dimension zs(2,2), zu(2,2), zui(2,2), zdef(2,2), zdel(2,2)
       dimension ztmp(2,2)
      common /kinemt/   El, Ecm, k0, pl, s, Mp2, Mnuc, Mnuc2
      common /params/   hbarc, pi, Mp, MN, nz, na, nes, nwaves
      common /ranges/   rcoul, rcut

            zi = (0.d0, 1.d0)
c      reset coul functions
            eta = 0.
            rhoc = k0*rcut/hbarc
            lcmax = ldmax + 2
            if (lcmax.le.1) lcmax = 2
            accur = 1.e-14
            step = 999
            nminl = 0
            call rcwfn (rhoc,eta,nminl,lcmax,fc0,fcp0,gc0,gcp0,accur,
     1      step)
c     correction of coulomb for a/(a-1)factor in kmt prescription
c           eta = xgam/aovera
            eta = xgam
            call rcwfn (rhoc,eta,nminl,lcmax,fc,fcp,gc,gcp,accur,step)
c          endif

c	jdum is actual total angular momentum starting from 0
 	if (ljtab .eq. 0) then
          ld = jdum 
          ld2 = ld + 2 
        else
          ld = jdum + 1
          ld2 = ld
        endif


c  *******  OK, now do the parameterization  *****************
            zt00 = dcmplx(tr00,ti00)
            zt02 = dcmplx(tr02,ti02)
            zt20 = dcmplx(tr20,ti20)
            zt22 = dcmplx(tr22,ti22)

c      write(7,1066) zt00, zt02
c      write(8,1066) zt22, zt20
c      write(8,1067) jdum, ljtab, zt22
c      write(7,1067) jdum, ljtab, zt00

c ************************************************************
            zs(1,1) = 1.d0 + 2.d0 * zi * zt00
            zs(1,2) = +2.d0 * zi * zt02
            zs(2,1) = +2.d0 * zi * zt20
            zs(2,2) = 1.d0 + 2.d0 * zi * zt22
            pud = dreal(zs(2,1))**2.d0 + dimag(zs(2,1))**2.d0
            if (pud.lt.1.e-16) write(8,*) "Problem 2"
            pud = dreal(zs(1,2))**2.d0 + dimag(zs(1,2))**2.d0
            if (pud.lt.1.e-16) write(8,*) "Problem 4"

            zrrr = zsqrt((zs(1,1) - zs(2,2))**2.d0 +
     &              4.d0 * zs(2,1) * zs(1,2))
            zlamp = 0.5 * (zs(1,1) + zs(2,2) + zrrr)
            zlamm = 0.5 * (zs(1,1) + zs(2,2) - zrrr)
            pud = dreal(zrrr)**2.d0 + dimag(zrrr)**2.d0
            if (pud.lt.1.e-16) write(8,*) "Problem 5"
            pud = dreal(zlamm)**2.d0 + dimag(zlamm)**2.d0
            if (pud.lt.1.e-16) write(8,*) "Problem 6"
            zdet = 1.d0 - (zs(1,2) * zs(2,1))
     &                    / ((zlamp - zs(1,1)) * (zlamm - zs(2,2)))
            zlet = 0.d0 - (zs(1,2) * zs(2,1))
     &                    / ((zlamp - zs(1,1)) * (zlamm - zs(2,2)))
            pud = dreal(zlet)**2.d0 + dimag(zlet)**2.d0
            if (pud.lt.1.e-16) write(8,*) "Problem 1"
            zud = zs(1,1) - zs(2,2)
            pud = dreal(zud)**2.d0 + dimag(zud)**2.d0
            if (pud.lt.1.e-16) write(8,*) "Problem 3"
c            write(7,1066) zdet, zlet
c           write(7,*) pi

            zu(1,1) = 1.d0
            zu(2,1) = zs(2,1) / (zlamm - zs(2,2))
            zu(1,2) = zs(1,2) / (zlamp - zs(1,1)) 
            zu(2,2) = 1.d0

            zui(1,1) = 1.d0/ zdet
            zui(1,2) = -(zs(1,2) / (zlamp - zs(1,1)))/zdet
            zui(2,1) = -(zs(2,1) / (zlamm - zs(2,2)))/zdet
            zui(2,2) = 1.d0/zdet

            do 111 i = 1, 2
            do 111 j = 1, 2
            zdef(i,j) = 0.0
            ztmp(i,j) = 0.0
            do 111 k = 1, 2
            ztmp(i,j) = ztmp(i,j) + zu(i,k) * zui(k,j)
            zdef(i,j) = zdef(i,j) + zui(i,k) * zs(k,j)
 111        continue
            zud = ztmp(2,1)
            pud = dreal(zud)**2.d0 + dimag(zud)**2.d0
            if (pud.gt.1.e-18) write(8,*) "Problem 10"


            do 222 i = 1, 2
            do 222 j = 1, 2
            zdel(i,j) = 0.0
            do 222 k = 1, 2
            zdel(i,j) = zdel(i,j) + zdef(i,k) * zu(k,j)
 222        continue
            pud = dreal(zdel(1,2))**2.d0 + dimag(zdel(1,2))**2.d0
            if (pud.gt.1.e-18) write(8,*) "Problem 7"


c ******** OK, now actually do the matching  *****************
             zd0 = zsqrt(zdel(1,1))
             zd2 = zsqrt(zdel(2,2))
             pud = dreal(zd0)**2.0 + dimag(zd0)**2.0 
             puddy = dreal(zd2)**2.0 + dimag(zd2)**2.0
c            if (pud.gt.1.0.or.puddy.gt.1.0) then
c            write(11,1139) ld, ljtab, pud, puddy
c            endif
c            if (pud.gt.1.0) then
c            puddy = dimag(zd0)
c            write(11,1138) ld, puddy
 1139        format(2i8,2e14.6)
 1138        format(i8,e14.6)
c            endif
c            if (pud.gt.e-0) write(8,*) "e^2idelta**2 > 1"
             ztan0 = -zi * (zd0 - (1.0 / zd0))/(zd0 + (1.0/zd0))
             ztan2 = -zi * (zd2 - (1.0 / zd2))/(zd2 + (1.0/zd2))
c          ******** Wronskians first ***************
       w1 = fcp(ld) * fc0(ld) - fcp0(ld) * fc(ld)
       w2 = fcp(ld) * gc0(ld) - gcp0(ld) * fc(ld)
       w3 = gcp0(ld) * gc(ld) - gcp(ld) * gc0(ld)
       w4 = fcp0(ld) * gc(ld) - gcp(ld) * fc0(ld)
       ww1 = fcp(ld2) * fc0(ld2) - fcp0(ld2) * fc(ld2)
       ww2 = fcp(ld2) * gc0(ld2) - gcp0(ld2) * fc(ld2)
       ww3 = gcp0(ld2) * gc(ld2) - gcp(ld2) * gc0(ld2)
       ww4 = fcp0(ld2) * gc(ld2) - gcp(ld2) * fc0(ld2)


        ztan0 = (w1 + ztan0 * w2) / (ztan0 * w3 + w4)
        ztan2 = (ww1 + ztan2 * ww2) / (ztan2 * ww3 + ww4)
            zud = ztan0 * w3 + w4
            pud = dreal(zud)**2.d0 + dimag(zud)**2.d0
            if (pud.lt.1.e-16) write(8,*) "Problem 8"

        zdel(1,1) = (1.0 + zi * ztan0)/(1.0 - zi * ztan0)
        zdel(2,2) = (1.0 + zi * ztan2)/(1.0 - zi * ztan2)
            zud = 1.0 - zi * ztan2
            pud = dreal(zud)**2.d0 + dimag(zud)**2.d0
            if (pud.lt.1.e-10) write(8,*) "Problem 9"
c       zdel(1,2) = 0.d0
c       zdel(2,1) = 0.d0

c ********OK, Now go back to the t-matrix    *****************
            do 333 i = 1, 2
            do 333 j = 1, 2
            zdef(i,j) = 0.0
            ztmp(i,j) = 0.0
            do 333 k = 1, 2
            zdef(i,j) = zdef(i,j) + zu(i,k) * zdel(k,j)
            ztmp(i,j) = ztmp(i,j) + zui(i,k) * zu(k,j)
 333        continue
            zud = ztmp(2,1)
            pud = dreal(zud)**2.d0 + dimag(zud)**2.d0
            if (pud.gt.1.e-18) write(8,*) "Problem 12"


            do 444 i = 1, 2
            do 444 j = 1, 2
            zs(i,j) = 0.0
            do 444 k = 1, 2
            zs(i,j) = zs(i,j) + zdef(i,k) * zui(k,j)
 444        continue


        zt00 = -zi / (2.d0) * (zs(1,1) -1.d0)
        zt02 = -zi * zs(1,2) / (2.d0) 
        zt20 = -zi * zs(2,1) / (2.d0)
        zt22 = -zi / (2.d0) * (zs(2,2) -1.d0)

c **********   print things out, see if they make sense ******
c      write(7,1066) zt00, zt02
c      write(8,1066) zt22, zt20
c      write(7,1066) zd0, zd2
c      write(8,1066) ztmp(2,2), ztmp(2,1)
c      write(8,1067) jdum, ljtab, zt22
c      write(7,1067) jdum, ljtab, zt00
 1066  format( 2('(',e16.8,',',e16.8,')'))
 1067  format(2i8, 1('(',e16.8,',',e16.8,')'))
 1068  format(i8, e16.8)
c      write(8,*) "                    "
c      write(7,*) "                    "

c ************************************************************

          if (ljtab.lt.0.5) then
c         write(7,1068) jdum, tr00
c         write(77,1068) jdum, ti00
          else
c         write(8,1068) jdum, tr00
c         write(88,1068) jdum,ti00
          endif

c
            tr00 = dreal(zt00)
            ti00 = dimag(zt00)
            tr02 = dreal(zt02)
            ti02 = dimag(zt02)
	    tr20 = dreal(zt20)
            ti20 = dimag(zt20)
            tr22 = dreal(zt22)
            ti22 = dimag(zt22)

          if (ljtab.lt.0.5) then
c         write(7,1068) jdum, tr00
c         write(77,1068) jdum, ti00
          else
c         write(8,1068) jdum, tr00
c         write(88,1068) jdum,ti00
          endif


	  return
      end

