c// Lu, 7/16/91, changed this routine for direct matching of t matrix.
c// Tim, 4/16/94, rewrote to match coupled channel case.
c// Input: 4 complex t matrix; output: same 4 complex t matrix
c// Blatt Biedenharn convention
c ***   See Lu, Mefford, Landau and Song, Phys. Rev. C 50, 3037, (1994)
      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 ************************************************************
            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

            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)
            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)))

            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

            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


c ******** OK, now actually do the matching  *****************
             zd0 = zsqrt(zdel(1,1))
             zd2 = zsqrt(zdel(2,2))
             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)

        zdel(1,1) = (1.0 + zi * ztan0)/(1.0 - zi * ztan0)
        zdel(2,2) = (1.0 + zi * ztan2)/(1.0 - zi * ztan2)

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

            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)

 1066  format( 2('(',e16.8,',',e16.8,')'))
 1067  format(2i8, 1('(',e16.8,',',e16.8,')'))
 1068  format(i8, e16.8)

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

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)


	  return
      end

