c$debug
c***********************************************************************
      subroutine ppbar(ii,jj,lk,ncup,icode,ippbar,x
     1                 ,space,gf,w,abeta0)
c
c     Subroutine to calculate coupling potentials for the scalar 
c     model of a p p_bar to hyperon hyperon_bar reaction
c
c	ii	initial channel
c	jj	final   channel
c	lk	number of present coupling matrix
c	ncup	number of present couplings
c	icode	form factor option number
c	ippbar	particle-hole code = true for p-h
c	x	array for form factor outputs
c	space	form factor shape storage (inputted)
c	gf	array for form factor calculation
c	w	optical potential storage
c	abeta0	strength parameter for form factor
c
c***********************************************************************
c
      parameter (kmax = 300, ispc0 = 16384, nlev = 12, nlsj = 9)
      parameter (nlkmax0 = ispc0/2/kmax, nlkmax = 3*nlkmax0)
c
      logical ippbar
      complex xtemp(nlsj), X(kmax,nlkmax0),W(kmax,8)
c
      COMMON ALPHA(15),HBARC,AMU,CHGSQ,DRT,RMAX
     1,ICON(20),LMAX,KT,MAXD,IOUT,JOUT,NLTR,IDAT(6)
     2,KLAM(nlev),LTRT(nlev),JTRT(nlev),NADD(nlev)
      COMMON EN(nlev),FM(nlev),Z(nlev),FMA(nlev),ZA(nlev),RY(nlev)
     1,FS(nlev),ECM(nlev),FMU(nlev),RC(nlev),FK2(nlev),FK(nlev)
     2,ETA(nlev),AFACT(nlev),PFACT(nlev),DRF(nlev),FN(nlev)
     3,PNLOC(nlev),DUM(nlev),TSIG(nlev)
     4,IS(nlev),NS(nlev),NLI(nlev),IN(nlev),INN(nlev,nlev)
     5,LAM(nlkmax),ISP(nlkmax),JSP(nlkmax),ILK(nlkmax)
     6,LKTAB(2*nlkmax),ICUP(2*nlkmax),JCUP(2*nlkmax)
	common/ch3files/input, ioutp
c
      dimension gf(2*kmax),space(8*kmax)
     1,c0(4),c1(3),c2(3),c4(4),c5(4),c6(4),lsj(3,nlsj)
     2,ilkl(16),iikl(16), h_mass(4), is12(2,3)
c
c     data for derivative coupling case
      data lsj/0,2,2
     2,        2,2,2
     3,        2,2,2
     4,        2,2,2
     5,        1,2,0
     6,        1,0,2
     7,        0,0,0
     8,        0,0,0
     9,        0,1,1/
c
c     masses for the proton, lambda_0, sigma_0 and sigma_+
      data h_mass/   1.0078, 1.1976,   1.2802,     1.2802/
c     passive quark spin couplings (lam,lam_bar), (lam,sig_bar),(sig,sig_bar)
      data is12/ 0, 0, 1,  0, 1, 1/
c
      i=iabs(ii)
      j=iabs(jj)
      abeta = abeta0
c
      if(ippbar) then
        ipp = 1
        if(abs(fm (i) - h_mass(1)) .le. 0.0010)then
c
c       initial state is a p p_bar, determine the type of reaction 
c
          if     (abs(fm (j) - h_mass(2)) .le. 0.0010 .and.
     1            abs(fma(j) - h_mass(2)) .le. 0.0010)then
c         final state is a lambda lambda_bar
            ks1 = is12(1,1)
            ks2 = is12(2,1)
          elseif (abs(fm (j) - h_mass(2)) .le. 0.0010 .and.
     1            abs(fma(j) - h_mass(3)) .le. 0.0010)then
c         final state is a lambda sigma_bar
            ks1 = is12(1,2)
            ks2 = is12(2,2)
          elseif (abs(fm (j) - h_mass(4)) .le. 0.0010 .and.
     1            abs(fma(j) - h_mass(4)) .le. 0.0010)then
c         final state is a sigma  sigma_bar
            ks1 = is12(1,3)
            ks2 = is12(2,3)
          else
            stop ' Error --- masses do not match known reaction'
          endif
c         a1 is ratio of strange quark mass to up quark mass
          a1 = 3.0*(fm (j)+fma(j))/(fm (i)+fma(i)) - 2.0
          a2 = -a1
          a4 = 1.0
          a3 = -a4
c         multiply by 1/(9*m*m')
          abeta = abeta*(hbarc/(amu*fm(i)))**2/a1
        else
          read (input,9008)a1,a2,a3,a4,fks1,fks2
          ks1 = fks1
          ks2 = fks2
        endif
c
      else
        ipp = 0
        read (input,9008)a1,a2,a3,a4,fks1,fks2
        ks1 = fks1
        ks2 = fks2
      endif
c
      if(icon(19).eq.0) then
        write(ioutp,9012)a1,a2,a3,a4
      endif
c         coefficients c_kSS'
	  c000  = coeff_fct(0,0,ks2,ks1,ipp,0)
	  c011  = coeff_fct(1,1,ks2,ks1,ipp,0)
          if(icode .eq. 21) then
	  c110  = coeff_fct(0,1,ks2,ks1,ipp,1)
	  c111  = coeff_fct(1,1,ks2,ks1,ipp,1)
	  c211  = coeff_fct(1,1,ks2,ks1,ipp,2)
          endif
c	write (0,'(a,8f12.6)')' c_kSS ',c000,c011,c110,c111,c211
c
         if(ippbar) then
          if(icon(19).eq.0) then
            write(ioutp,9013)
          endif
          c000 = -2.0*c000
          c011 = -2.0*c011
        endif
c
      IF(icode .eq. 20) THEN
      c0(1) = abeta*(0.25*c011+0.75*c000)
      c6(1) = abeta*(0.75*c011-0.75*c000)
      do 2020 m=1,kmax
      mk = 2*m-1
c     L=0, sig*sig term
      x(m ,lk  ) = c0(1)*space(mk)
c     L=0, scalar
      x(m ,lk+1) = c6(1)*space(mk)
 2020 continue 
      ilk(lk) = 8
      if(icon(19).eq.0) then
        write(ioutp,9037)lam(lk),isp(lk),jsp(lk),ilk(lk)
      endif
      lk=lk+1
      ilk(lk) = 8
      ncup = ncup+1
      icup(ncup) = i
      jcup(ncup) = j
      lktab(ncup) = lk
      inn(j,i) = 1
      if(j.ne.i .and. jj.gt.0) then
        ncup = ncup+1
        icup(ncup) = j
        jcup(ncup) = i
        lktab(ncup) = lk
        inn(i,j) = 1
      endif
      lam(lk) = 0
      isp(lk) = 0
      jsp(lk) = 0
      if(icon(19).eq.0) then
        write(ioutp,9037)lam(lk),isp(lk),jsp(lk),ilk(lk)
      endif
c     option 21
      ELSE
      Recoil = (drf(j)/drf(i))**2
      c0(1)= (abeta*(a3+a1-a2)/6.0        )*(0.25*c011+0.75*c000)
      c0(2)= (abeta*(a2+a1-a3)*Recoil/6.0 )*(0.25*c011+0.75*c000)
      c0(3)= (abeta*(a1-a2-a3+2*a4)/6.0   )*(0.25*c011+0.75*c000)
      c0(4)=-(abeta* a4                   )*(0.25*c011+0.75*c000)
c	write (0,'(a,8f12.6)')' c0(i) ',c0,Recoil
      c1(1)= (abeta*(a3+a1-a2)/6.0        )* c211
      c1(2)= (abeta*(a2+a1-a3)*Recoil/6.0 )* c211
      c1(3)= (abeta*(a1-a2-a3+2*a4)/6.0   )* c211
      c2(1)=-(abeta*a1/12.0               )* c211
      c2(2)= (abeta*a1*Recoil/12.0        )* c211
      c2(3)= (abeta*(a3-a2)/12.0          )* c211
      c31  = (abeta*a1/12.0               )* c211
      c4(1)= (abeta*sqrt(3.0)*(a1-a2-a3)  )*(c111+c110)
      c5(1)= (abeta*sqrt(3.0)*(a1-a2-a3)  )*(c111-c110)
      c6(1)= (abeta*(a3+a1-a2)/6.0        )*(0.75*c011-0.75*c000)
      c6(2)= (abeta*(a2+a1-a3)*Recoil/6.0 )*(0.75*c011-0.75*c000)
      c6(3)= (abeta*(a1-a2-a3+2*a4)/6.0   )*(0.75*c011-0.75*c000)
      c6(4)=-(abeta* a4                   )*(0.75*c011-0.75*c000)
c
      r1 =drf(i)*float(kmax+1)
      t2 =0.0
      gft=0.0
      do 2030 m=2,kmax
      mm=   kmax+2 -m
      mk=(2*kmax+3)-2*m
      r2=r1
      r1=r2-drf(i)
c     gf is integral of H/R
c     h2 is 2nd derivative of H
c     h1 is 1st derivative of H
      gft=gft-(t2/r2+space(mk)/r1)*(drf(i)*0.5)
      gf(mk)=gft
      h2 =(t2-2.0*space(mk)+space(mk-2))/drf(i)**2
      h1 =(t2-space(mk-2))/(2.0*drf(i))
      t2=space(mk)
c	if(mod(mk+1,10).eq.0) write(21,7777)r1,space(mk),h1,h2,gft
 7777	format(f8.2,1p10e12.4)
c
c     L=0, sig*sig term
      x(mm ,lk  )=t2*(c0(1)*w(mm ,i) + c0(2)*w(mm ,j))
     1   + c0(3)*(h2+2.0*h1/r1)
c      if(correct.ne.0.0) then
c        x(mm ,lk  )=x(mm ,lk  ) + c0(4)*h1/r1
c      endif
c
c	if(mm.eq.5) write(21,7777)r1,space(mk),gft
c     1	,c0(1),w(mm,i),c0(2),w(mm,j),x(mm,lk)
c
c     L=2, tensor term
      x(mm ,lk+1)=t2*(c1(1)*w(mm ,i) + c1(2)*w(mm ,j))
     1   + c1(3)*(h2-    h1/r1)
c
c     tensor * (l**2)    operator
      x(mm ,lk+2)=gft*(c2(1)*w(mm ,i) + c2(2)*w(mm ,j))
     1   + c2(3)*(h1-2.*t2/r1)/r1
c
c     tensor * (l**2)**2 operator
      x(mm ,lk+3)=    c31*(gft-0.5*t2)/r1**2
c
c     spin orbit term, projectile
      x(mm ,lk+4)=    c4(1)*h1/r1
c
c     spin orbit term, target
      x(mm ,lk+5)=    c5(1)*h1/r1
c
c     L=0, scalar
      x(mm  ,lk+6)=t2*(c6(1)*w(mm ,i) + c6(2)*w(mm ,j))
     1   + c6(3)*(h2+2.0*h1/r1)
c      if(correct.ne.0.0) then
c        x(mm ,lk+6)=x(mm ,lk+6) + c6(4)*h1/r1
c      endif
c
c     scalar  * (l**2)    operator
      x(mm ,lk+7)=    0.0
c
c     sig*sig * (l**2)    operator
      x(mm ,lk+8)=    0.0
c
 2030 continue
c
c     extrapolate to obtain first point
      do 2035 k=1,nlsj
      x(1,lk+k-1)= (8.0*x(2,lk+k-1)-3.0*x(3,lk+k-1))/5.0
 2035 continue
      gf(1)=(8.0*gf(3)-3.0*gf(5))/5.0
c
      lk=lk-1
      ncup=ncup-1
      if(j.ne.i.and.jj.gt.0) then
        ncup=ncup-1
      endif
      do 2040 k=1,nlsj
      lk=lk+1
      ilk(lk)=icode+k-1
      ncup=ncup+1
      icup(ncup)=i
      jcup(ncup)=j
      lktab(ncup)=lk
      inn(j,i)=1
      if(j.ne.i.and.jj.gt.0) then
        ncup=ncup+1
        icup(ncup)=j
        jcup(ncup)=i
        lktab(ncup)=lk
        inn(i,j)=1
c
      endif
      lam(lk)=LSJ(1,K)
      isp(lk)=LSJ(2,K)
      jsp(lk)=LSJ(3,K)
      if(icon(19).eq.0) then
        write(ioutp,9037)lam(lk),isp(lk),jsp(lk),ilk(lk)
      endif
 2040 continue
c
c     contributions from non-scalar diagonal potentials
c
      lk1 = lk-nlsj
      if(lk1.gt.0) then
      nlkl=0
      do 2045 k=1,lk1
      if(ilk(k).eq.7.and.((icup(k).eq.I.and.jcup(k).eq.I)
     1 .or.(icup(k).eq.J.and.jcup(k).eq.J))) then
        nlkl=nlkl+1
        ilkl(nlkl)=k
        if(icup(k).eq.I) then
          iikl(nlkl)=1
        else
          iikl(nlkl)=2
        endif
      endif
 2045 continue
c
      if(nlkl.ne.0) then
      do 2090 m=1,kmax
      r1 = float(m)*drf(i)
      mk = 2*m-1
      do 2050 n=1,nlsj
      xtemp(n  )=0.0
 2050 continue
      do 2070 n=1,nlkl
      k =ilkl(n)
      i =icup(n)
      i1=iikl(n)
        if    (lam(k).eq.0.and.isp(k).eq.0) then
c
c     Scalar
c
          xtemp( 1)=xtemp( 1)+    c0(i1)*x(m  ,k)*space(mk)*fmu(i)
          xtemp( 2)=xtemp( 2)+    c1(i1)*x(m  ,k)*space(mk)*fmu(i)
          xtemp( 3)=xtemp( 3)+    c2(i1)*x(m  ,k)*gf   (mk)*fmu(i)
          xtemp( 7)=xtemp( 7)+    c6(i1)*x(m  ,k)*space(mk)*fmu(i)
c
c	if(m.eq.5) write(21,7777)r1,space(mk),gf(mk)
c     1	,c0(i1),x(m,k),x(m,k)*fmu(i)*gf(mk),xtemp(1)
c
        elseif(lam(k).eq.0.and.isp(k).eq.2) then
c
c     Sigma-Sigma
c
          cx=c6(i1)-2.0*c0(i1)
          xtemp( 1)=xtemp( 1)+    cx    *x(m  ,k)*space(mk)*fmu(i)
          xtemp( 2)=xtemp( 2)+    c1(i1)*x(m  ,k)*space(mk)*fmu(i)
          xtemp( 3)=xtemp( 3)+    c2(i1)*x(m  ,k)*gf   (mk)*fmu(i)
          xtemp( 7)=xtemp( 7)+3.0*c0(i1)*x(m  ,k)*space(mk)*fmu(i)
c
        elseif(lam(k).eq.2.and.isp(k).eq.2) then
c
c     Tensor
c
          cx=c0(i1)+c6(i1)-2.0*c1(i1)
          xtemp( 2)=xtemp( 2)+    cx    *x(m  ,k)*gf   (mk)*fmu(i)
          xtemp( 3)=xtemp( 3)-2.0*c2(i1)*x(m  ,k)*gf   (mk)*fmu(i)
c
          if(ippbar) then
            xtemp( 1)=xtemp( 1)+2.0*c1(i1)*x(m  ,k)*space(mk)*fmu(i)
            xtemp( 7)=xtemp( 7)+6.0*c1(i1)*x(m  ,k)*space(mk)*fmu(i)
            xtemp( 8)=xtemp( 8)+6.0*c2(i1)*x(m  ,k)*gf   (mk)*fmu(i)
            xtemp( 9)=xtemp( 9)+2.0*c2(i1)*x(m  ,k)*gf   (mk)*fmu(i)
          else
            xtemp( 7)=xtemp( 7)+8.0*c1(i1)*x(m  ,k)*space(mk)*fmu(i)
            xtemp( 8)=xtemp( 8)+8.0*c2(i1)*x(m  ,k)*gf   (mk)*fmu(i)
          endif
c
        endif
c	if(m.eq.5) write(21,7778)n,k,i1,lam(k),isp(k),space(mk),x(m  ,k)
c     1	,xtemp(1),xtemp(2),xtemp(3)
 7778	format(5i4,1p10e12.4)
 2070   continue
        do 2080 n=1,nlsj
        x(m  ,lk-nlsj+n)=x(m  ,lk-nlsj+n)+xtemp( n)
 2080   continue
c
 2090 continue
      endif
      endif
      ENDIF
c
      return
c
 9008 FORMAT(10F8.4)
 9012 format('0Coefficients        A1   =',f9.4,'   A2   =',f9.4
     1,       '   A3   =',f9.4,'   A4   =',f9.4 /)
 9013 format('0Particle-hole projection enabled')
 9037 FORMAT(18H                  ,9H   LAM  =,I4,14H        2*STR=,I4
     1,14H        2*JTR=,I4,14H        ICODE=,I4)
      end

c***********************************************************************
      subroutine gamma(ii,jj,lk,ncup,icode,ippbar,x,space,gf,w,abeta)
c
c     Subroutine to calculate coupling potentials for a (p,gamma) or
c     (gamma,p) reaction.   ***Not finished***
c     
c***********************************************************************
c
      parameter (kmax = 300, nlsj= 9, nlkmax0 = 16384/2/kmax, nlev = 12)
      parameter (nlkmax = 2*nlkmax0)
      logical ippbar
c
      COMMON ALPHA(15),HBARC,AMU,CHGSQ,DRT,RMAX
     1,ICON(20),LMAX,KT,MAXD,IOUT,JOUT,NLTR,IDAT(6)
     2,KLAM(nlev),LTRT(nlev),JTRT(nlev),NADD(nlev)
      COMMON EN(nlev),FM(nlev),Z(nlev),FMA(nlev),ZA(nlev),RY(nlev)
     1,FS(nlev),ECM(nlev),FMU(nlev),RC(nlev),FK2(nlev),FK(nlev)
     2,ETA(nlev),AFACT(nlev),PFACT(nlev),DRF(nlev),FN(nlev)
     3,PNLOC(nlev),DUM(nlev),TSIG(nlev)
     4,IS(nlev),NS(nlev),NLI(nlev),IN(nlev),INN(nlev,nlev)
     5,LAM(nlkmax),ISP(nlkmax),JSP(nlkmax),ILK(nlkmax)
     6,LKTAB(2*nlkmax),ICUP(2*nlkmax),JCUP(2*nlkmax)
	common/ch3files/input, ioutp
c
      dimension X(2*kmax,32),W(2*kmax,8),space(8*kmax),gf(2*kmax)
c
      if(ippbar) i=0
      i=iabs(ii)
      j=iabs(jj)
c
      c1 = abeta
      r1 = drf(i)*float(kmax+1)
      r0 = drf(i)*float(kmax)
      t2 =0.0
      gft=0.0
      do 1620 m=2,kmax
      mk=(2*kmax+3)-2*m
      r2=r1
      r1=r0
      r0=r1-drf(i)
c     gf is integral of R*H
c     h1 is 1st derivative of R*H
      gft=gft-(t2*r2+space(mk)*r1)*(drf(i)*0.5)
      gf(mk)=gft
      h1 =(t2*r2-space(mk-2)*r0)/(2.0*drf(i))
      t2=space(mk)
c	if(mod(mk+1,20).eq.0) write(21,7777)r1,space(mk),h1,gft
 7777	format(f8.2,1p10e12.4)
c     Scalar term
      x(mk  ,lk  )=gft*c1*(w(mk  ,i) - w(mk  ,j))
      x(mk+1,lk  )=gft*c1*(w(mk+1,i) - w(mk+1,j))

c     Scalar * (l**2)    operator
      x(mk  ,lk+2)=h1*c1/r1**2
      x(mk+1,lk+2)=0.0
 1620 continue
c     extrapolate to obtain first point
      do 1630 k=1,2
      x(1,lk+k-1)= (8.0*x(3,lk+k-1)-3.0*x(5,lk+k-1))/5.0
      x(2,lk+k-1)= (8.0*x(4,lk+k-1)-3.0*x(6,lk+k-1))/5.0
 1630 continue
      gf(1)=(8.0*gf(3)-3.0*gf(5))/5.0
c
      lk1=lk
      lk=lk-1
      ncup=ncup-1
      IF(J.NE.I.AND.JJ.GT.0) then
        ncup=ncup-1
      endif
      DO 1640 k=1,2
      lk=lk+1
      ilk(lk)=icode+k-1
      ncup=ncup+1
      icup(ncup)=i
      jcup(ncup)=j
      lktab(ncup)=lk
      inn(j,i)=1
      if(j.ne.i.and.jj.gt.0) then
        ncup=ncup+1
        icup(ncup)=j
        jcup(ncup)=i
        lktab(ncup)=lk
        inn(i,j)=1
c
      endif
      LAM(LK)=lam(lk1)
      ISP(LK)=isp(lk1)
      JSP(LK)=jsp(lk1)
      if(icon(19).eq.0) then
        WRITE(ioutp,9037)LAM(LK),ISP(LK),JSP(LK),ILK(LK)
      endif
 1640 continue
c
      return
c
 9037 FORMAT(18H                  ,9H   LAM  =,I4,14H        2*STR=,I4
     1,14H        2*JTR=,I4,14H        ICODE=,I4)
      end

c$debug
	function coeff_fct(jS,iS,js12,is12,ipp,k)
c
c	calculates factors for conversion from microscopic q=q_bar
c	to macroscopic hyperon-hyperon_bar interaction
c
c	jS	final   state composite particle spin
c	iS	initial state composite particle spin
c	js12	second  particle spectator quark spin
c	is12	first   particle spectator quark spin
c	kS12	spin coupling of spectator quarks
c	ipp	= 0,  particle-particle case
c		= 1,  particle-hole case
c	k	rank of interaction spin tensor
c
c	kS12	spin coupling of spectator quarks
c	js3	final   state active quark spin
c	is3	initial state active quark spin
c
	coeff = 0.0
	if(phasef(iS  +k -jS ) .gt. 0.0 .or. ipp .eq. 1)then
	  k1 = 1
	  k2 = 1
	if(ipp .eq. 1 .and. k .eq. 1)then
	  k2 = 0
	endif
	n1 = abs(is12-js12)
	n2 =     is12+js12
	do 500 kS12= n1 ,n2
	i1 = max(abs(iS-kS12), ipp)
	i2 = min(    iS+kS12,  1  )
	do 400 is3 = i1, i2
	j1 = max(abs(jS-kS12), ipp)
	j2 = min(    jS+kS12,  1  )
	do 300 js3 = j1, j2
	if(k    .ge. abs(is3-js3) .and. k    .le. is3+ js3) then
c
	Cjj_LS = sqrt(float((2*kS12+1)*(2*is3+1)*4))
     1         * winej(2*is12, 2*js12, 2*kS12, 1, 1, 2*is3, 1, 1, 2*iS)
     2         * sqrt(float((2*kS12+1)*(2*js3+1)*4))
     3         * winej(2*is12, 2*js12, 2*kS12, 1, 1, 2*js3, 1, 1, 2*jS)
c
	RME1   = sqrt(float((2*is3+1)*(2*js3+1)*(2*k+1)))*6.0
     1         * winej(1, 1, 2*k1, 1, 1, 2*k2, 2*js3, 2*is3, 2*k)
	RME2   = sqrt(float((2*iS +1)*(2*jS +1)*(2*k+1)))*6.0
     1         * winej(1, 1, 2*k1, 1, 1, 2*k2, 2*jS , 2*iS , 2*k)
c
	coefft = phasef(k+kS12-iS-js3)
     1	       * racah(2*is3,2*iS,2*js3,2*jS,2*kS12,2*k)
     2	       * sqrt(float((2*iS +1)*(2*jS +1)))
	coeff = coeff + coefft *Cjj_LS *RME1 /RME2 
c	write(6,7777)Cjj_LS,RME1,RME2,coefft
c 7777	format(' Cjj_LS, RME1, RME2, coefft = '  ,4f12.4)
	endif
  300	continue
  400	continue
  500	continue
	endif
	coeff_fct = coeff
	return
	end
