c/* @(#)tmatrx.f	1.1      latest revision 6/23/89 14:55:21 */
      subroutine coulomb( scoul, sigl, aovera, ld, nspin, njump )
**********************************************************************

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

      complex*16  zi , ztan
      complex*16  zr34 , zr56 , zrmix , ztjp, ztjm , ztj

      real*8      k0
      real*8      Mnuc , Mnuc2 , MN , Mp , Mp2

c *** N.B. need to match dimensions and # of entries in data statement
c *** for the IBM.
c *** t's enlarged 3/84 for cex  (m.s.)

      dimension tbr(100,12), tbi(100,12), tr(100,12),ti(100,12)

c *** dimension of u/f depends on ngp (only need change in main)
c *** double size for spin 0 case,so can up n grid pts

      dimension   scoul(100)
      dimension   nifty(20)

      common /kinemt/   El , Ecm , k0, pl , s , Mp2 , Mnuc , Mnuc2
      common /params/   hbarc, pi, Mp, MN, nz, na, nes, nwaves
      common /switch/   nifty
      common /sizes/    achp, acmp, wsp, achn, acmn, wsn
      common /ranges/   rcoul, rcut
      common /inputs/   Tlab , b , ymin1 , ymin2 , xnang,
     $                  kode , lxmax , ngp , nr

      common /Rcomn/    rr , ri , rrb , rib
      common /spins/    nspina , nspinb , nspinc , nspind , nspine
      common /Tcoul/    trc , tic, tr3c, ti3c, tr5c, ti5c, xgam
      common /Tcomn/    tr , ti , tbr , tbi

***************************************************************************

c *** calculation  of    t-matrix  (on - shell),
c *** t normalized to exp(i del)*sin(del)

      n = nspin
      zi = ( 0. , 1. )
      if ( njump .eq. 1 ) go to 290
      if ( nspin .ne. 4) go to 281
      go to 281

c ***
c *** stapp convention with blatt-beidenharn phases
c *** rhl change to nspin coding,l-2 subtractn in thepp
c *** ztjp=alpha(j+1,j)=t+     ztjm=alpha(j-1,j)=t-   ztj= alpha**(j+1)
c ***
c *** rhl version,reordered see nspin code and thep sub for stapp convrson

c     zr34 = rr + zi * ri
      go to 300
 281  if (nspin .ne. 6) go to 282
c     zr56 = rr + ri * zi
c     zrmix = rrb + zi * rib

c *** call rtot ,t(4) and t(6) now set =, the diffence is insig

      iborn = 0
      if (nifty(3)   .eq.  1
     $               .OR.
     $    nifty(3)   .eq.  2
     $               .OR.
     $    nifty(6)   .eq.  5)    iborn = 1

      write(6,3335) (tr(ld,l9), ti(ld,l9), l9=3,6)
 282  continue


      trc = tr(ld,n-1)
      tic = ti(ld,n-1)

c let run through phase mods with zero sigl
      if (nifty(10) .eq. 0) go to 300
      if (nifty(10) .lt. 3) go to 290
c if nifty(10) eq 3 or 4 the do exact coulomb (either uniform sphere (=3)
c or realistic (=4)


c *** ---- --coulomb v included on poten
c *** determine phase shift by matching
c ***
c ***    determine tan delta(sl) from calculated t*s

c     change 9/90 ldmax = lxmax
      ldmax = lxmax
******* tvm 10/22/92 change to all calls just nspin-1 **************
*********** Halloween only tvm 10/30/92 change *******************
c        write(11,666) ld, n-1, tr(ld,n-1), ti(ld,n-1)
c666     format(2i4, 2e14.7)
         jt = ld - 1
         jd = ld
***************** check this out ***********************************
         if(n.eq.8) then
         if(nifty(4).eq.1) then
         sssmall = 1.e-10
         tr(ld,8) = sssmall
         ti(ld,8) = sssmall * 1.0002
         tr(ld,2) = sssmall * .999
         ti(ld,2) = sssmall * 1.001
         endif
         if(ld.lt.3) then
         if(ld.eq.1) then
c               (ldum = 1 case here)
         singr = tr(1,1)
         singi = ti(1,1)
         ljtab = 0
         call matchc2(tr(ld,3),ti(ld,3),tr(ld,6),ti(ld,6),tr(ld,4),
     $     ti(ld,4),tr(ld,5),ti(ld,5),jd,ldmax,aovera,xgam,ljtab)
         call matchc(singr,singi,tr(1,1),ti(1,1),ld,ldmax,aovera,
     $              xgam,scoul)
         else
c             (ldum = 2 case here)   
         ljtab = 1
         call matchc2(tr(ld,1),ti(ld,1),tr(ld,2),ti(ld,2),tr(ld,8),
     $      ti(ld,8),tr(ld,7),ti(ld,7),jt,ldmax,aovera,xgam,ljtab)
          tugr = tr(1,7)
          tugi = ti(1,7)
         call matchc(tugr,tugi,tr(1,7),ti(1,7),ld,ldmax,aovera,
     $              xgam,scoul)
         ljtab = 0
         call matchc2(tr(ld,3),ti(ld,3),tr(ld,6),ti(ld,6),tr(ld,4),
     $     ti(ld,4),tr(ld,5),ti(ld,5),jd,ldmax,aovera,xgam,ljtab)
         endif
         else
c            (General case here)
         ljtab = 1
         call matchc2(tr(ld,1),ti(ld,1),tr(ld,2),ti(ld,2),tr(ld,8),
     $      ti(ld,8),tr(ld,7),ti(ld,7),jt,ldmax,aovera,xgam,ljtab)
         ljtab = 0 
         call matchc2(tr(ld,3),ti(ld,3),tr(ld,6),ti(ld,6),tr(ld,4),
     $     ti(ld,4),tr(ld,5),ti(ld,5),jd,ldmax,aovera,xgam,ljtab)
         endif
         endif
************* ugly, huh? ******************************************
c        write(11,666) ld, n-1, trc, tic
 290  continue
      njump = 0

c *** include coulomb phase factor in nuclear amplitude

      if(n.eq.8) then
      sigl = scoul(ld)
      if (nifty(10) .eq. 2) sigl = 0
      rhl5 = cos(2 * sigl)
      rhl6 = sin(2 * sigl)
      tr(ld,1) = rhl5 * tr(ld,1) - rhl6 * ti(ld,1)
      ti(ld,1) = rhl6 * tr(ld,1) + rhl5 * ti(ld,1)
      tr(ld,3) = rhl5 * tr(ld,3) - rhl6 * ti(ld,3)
      ti(ld,3) = rhl6 * tr(ld,3) + rhl5 * ti(ld,3)
       if (ld.gt.2) then
      tr(ld-2,5) = rhl5 * tr(ld-2,5) - rhl6 * ti(ld-2,5)
      ti(ld-2,5) = rhl6 * tr(ld-2,5) + rhl5 * ti(ld-2,5)
        endif
        if(ld.gt.1) then
      tr(ld,7) = rhl5 * tr(ld,7) - rhl6 * ti(ld,7)
      ti(ld,7) = rhl6 * tr(ld,7) + rhl5 * ti(ld,7)
       if (ld.lt.3) then
      tr(ld-1,7) = rhl5 * tr(ld-1,7) - rhl6 * ti(ld-1,7)
      ti(ld-1,7) = rhl6 * tr(ld-1,7) + rhl5 * ti(ld-1,7)
      endif
      endif
      endif

 300   continue

      RETURN

C *** formats

 1010 format (1h+,65x,4heta=,f10.4,8h  delta=,f10.4,'  sigl=',e14.6)
 1020 format (41h0********** eta greater than 1 **********,//)
 1030 format (1h ,15h  t(k,k,k)=    ,35x,4e14.7)
 1050 format (20x,40hcoulomb phase modified nuclear amplitude/16h   t(k,
     $k.k) =   ,4e14.7)
 3335 format(1h ,19h tr ti (nspin=3-6)=,8e13.5)

      end