Subroutine vabs (k,kp,k0,revl,imvl,nls)
c *** see RH Landau, Program lpott, 1981
c     sub to calc partial wave decomp of absortive contrib to opt pot
c       nifty(7)=3, the old case, i.e.
c     use spuller-measday params for pid-pp for e-dependence
c       now extented(6-77) to go beyond resonance
c     use empirical(atomic)values for zero energy (Read in via /absp/)
c     nifty(7)=4        use isobar model with different rhon,rhop possib
c     for he3 h3 special treatment with spin etc possible
c     then,first 3 params in sec4 are for magnetic, and last 3 for rho2
c     nb these  are a(mccarthy def)-p,b(mc carth def)-p, a(mcc)-n
c       the a s differ by factopr of 2 from  h.o. param
c     , for given value of k0,the params are fixed
      Implicit Real*8 (a-h,i,k,m,o-z)
      Real*8 xcos(48),wt(48),revl(30,2),imvl(30,2),pl(30)
      Dimension nifty(20), bnuc(20), bnucf(20), bn(16), bnf(16), retij
     1(14), imtij(14), ff(4), tofth(8), revofq(2), imvofq(2)
      Common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,mn,nz,nes
     1,nwaves,nifty,na
      Common /sec4/ achp,acmp,wsp,amccp,bmccp,amccn/absp/b0r,b0i,c0r,c0i
c
c     statement functions
c
      data kold,md/0.,1875.63/
      kcm2(ss) = (ss-(mpi+mn)**2)*(ss-(mpi-mn)**2)/4./ss
      eint(p1,p2) = Sqrt(mn2+(rhl1**2)*p1*p1*0.25+(rhl2**2)*p2*p2*0.25-0
     1.5*rhl1*rhl2*p1*p2*x)
      softh(p1,p2) = mpi2+mn2+rhl1*p1*p1-rhl2*p1*p2*x+2.*Sqrt(mpi2+p1*p1
     1)*eint(p1,p2)
      gs(p) = 776./(42.e04+p*p)
      gp(p) = 6.6e04*p/(8.18e04+p*p)**2
      rho2f(a2,alpha) = exp(-q2*a2/8.)*(1.+alpha*a2*(12.-q2*a2)/8.+(
     1(alpha*a2)**2)*(15.-2.5*q2*a2+0.0625*(q2*a2)**2)/16.)
      If (abs(1.-kold/k0).lt.1.e-03) GoTo 80
c     first cal,set up constants and pts wts
c     absap only for pions
      If (mpi.gt.140.) STOP
      rhl1 = (na+1.)/na
      rhl2 = (na-1.)/na
      rhl = na*(na-1.)
      rhl5 = na*na
      mpi2 = mpi*mpi
      mn2 = mn*mn
      Do 10 n=1,30
         Do 10 j=1,2
            revl(n,j) = 0.
            imvl(n,j) = 0.
 10   Continue
      ppth = Sqrt(0.25*(mpi2+md*md)+0.5*md*mpi-mn2)/mpi
c     set up pts for costheta projection
      ncoss = 16
      kold = k0
      Call Gauss (ncoss,011,-1.d0,1.d0,xcos,wt)
      b = .6*(1.1/hbarc)*(na**(.33333))
      If (nifty(16).ne.2) b = (achp+amccp)/2./hbarc
c       the rho**2 param are p and n average
      If ((na.eq.4).or.(na.eq.3)) b = 1.362/hbarc
      If (na.ge.100) b = 3./hbarc
      b2 = b*b
      alph = (na-4.)/6./b2
      If (alph.lt.0.) alph = 0.
c     normalization for rho2 when rho(r) nornmalized to one
      anrho2 = 0.5/(pi**(1.5))/(b2*b)/(1.+1.5*alph*b2)**2/Sqrt(2.)
      epi = Sqrt(mpi2+k0*k0)
      enuc = Sqrt(mn*mn*na*na+k0*k0)
c     ftot If the reduced mass factor,need pi2 in lpt norm
c     the number of pairs is now seperate and mult latter
      ftot = epi*enuc/(epi+enuc)
      anlpt = -4.*pi*pi*pi/2./ftot/(2.*pi)**3
      If (na.eq.13) rhl = rhl*12./13.
c     special calc for wood-saxon,Do rho2 numerically
      If (nifty(16).eq.2) anrho2 = 1.
      anorm = anrho2*anlpt*rhl
      jmax = 1
      If (nifty(6).eq.3) jmax = 2
c     block to calc energy dependce of b0 c0
c
c     calc pi d com momentum kpid,and p p com momentum pp
      If (nifty(7).eq.4) GoTo 70
      ksave = k0
      If (k0.lt.200.) GoTo 20
c     ko gt low e expansion range,eval at kcut then scale
      k0 = 200.
 20   ko = k0
      ed = Sqrt(md*md+(2.*k0/na)**2)
      epi = Sqrt(mpi2+k0*ko)
      ftot = epi*ed/(epi+ed)
      s = (epi+ed)**2-(k0*(1.-2./na))**2
      kpid2 = (s-(mpi+md)**2)*(s-(mpi-md)**2)/4./s/mpi2
      kpid = Sqrt(kpid2)
      pp = Sqrt(0.25*s-mn2)/mpi
c     s wave energy depepndce
      ftotp = mpi*md/(mpi+md)/ftot
      imb = ((b0i)/(mpi**4))*((pp/ppth)**2)*ftotp
      imb = imb*(1.-.81*kpid+0.31*kpid2*(1.-kpid2)+0.505*kpid**3)
      reb = b0r*imb/(b0i+1.e-08)
c     pwave,,, nr=non resonant contribution
      imc = (c0i/mpi**6)*ftotp*((pp/ppth)**2)*(1.+1.67*kpid-kpid2)
      rec = c0r*imc/(c0i+1.e-08)
      If (ksave.lt.200.) GoTo 60
c      k0 gt kcut, use reesonace form for re im p,polynml for s
      Do 50 n=1,2
         If (n.eq.1) GoTo 30
c     2nd time thru, redefine and scale
         k0 = ksave
         ko = k0
         ed = Sqrt(md*md+(2.*k0/na)**2)
         epi = Sqrt(mpi2+k0*k0)
         ftot = epi*ed/(epi+ed)
         s = (epi+ed)**2-(k0*(1.-2./na))**2
         kpid2 = (s-(mpi+md)**2)*(s-(mpi-md)**2)/4./s/mpi2
         kpid = Sqrt(kpid2)
         pp = Sqrt(0.25*s-mn2)/mpi
c     first time thru or second
 30      gamel = 0.6*((kpid)**3/(1.+kpid2))
         gamr = 0.05*0.92*pp
         gamt = 0.92*(gamel/0.6+0.05*pp)
         rhl5 = (5./3.)*pi*((mpi*kpid)**(-2))*gamel*gamr/(((Sqrt(s)-2179
     1   .)/mpi)**2+gamt*gamt/4.)
         rhl5 = rhl5*10.*hbarc**2
c     scale s wave up
         rhl5 = rhl5*4
         rhl6 = kpid*rhl5/4.+(2./3.)*pp*pp*(.26-.4*kpid+.08*kpid2)
         If (rhl6.lt.0.) rhl6 = 0.
         If ((k.eq.kp).and.(k.eq.k0)) Write (6,200) n,rhl5,rhl6
         If (n.eq.2) GoTo 40
c     1st time thru
         imc = imc/(rhl5*ftot/kpid)
         imb = imb/(rhl6*ftot)
         GoTo 50
c     2nd tyme thru
 40      imc = imc*ftot*rhl5/kpid
         imb = imb*ftot*rhl6
         reb = imb*b0r/(b0i+1.e-08)
         rec = imc*c0r/(c0i+1.e-08)
 50   Continue
 60   Continue
      GoTo 80
c     isobar calculation
c       can put partial wave decomp here,instead us
c     it as a function of theta,ie call tnuc for angle tf
c     intialize
c     no pairs for absp here, related to def of b0,c0,nb 1/4
 70   nn = na-nz
      pairpp = 0.
      pairnn = 0.
      pairpn = nz*nn
      If (nifty(1).eq.1) pairpp = nz*(nz-1)/2.
      If (nifty(1).eq.0) pairnn = nn*(nn-1)/2.
      If (na.ne.3) GoTo 80
c     special he3/h3 density part ( b refers to ho param)
      bp = amccp*2
      bp2 = bp*bp
      bmccp2 = bmccp*bmccp/bp2
      alphp = (bmccp2*4.)/(1.-6.*bmccp2)
      arhopp = 0.5/(pi**(1.5))/(bp2*bp)/(1.+1.5*alphp*bp2)**2/Sqrt(2.)
      alphn = 0.
      bbn = amccn*2.
      bbn2 = bbn*bbn
      bpn2 = bbn2*bp2/(bbn2+bp2)
      arhonn = 0.5/(pi**(1.5))/(bbn2*bbn)/Sqrt(2.)
      arhopn = 1./(bbn2*bbn)/(1.+1.5*alphp*bp2)/(2.*pi*(bp2+bbn2)/2./
     1bbn2)**1.5
c     put lpt norm and no pairs and rho2 norm togerther
      arhopp = arhopp*pairpp*anlpt
      arhonn = arhonn*pairnn*anlpt
      arhopn = arhopn*pairpn*anlpt
c     params for vabs all stored for this k0, now calc vabs
c
 80   Do 90 n=1,nls
         Do 90 j=1,jmax
            revl(n,j) = 0.
            imvl(n,j) = 0.
 90   Continue
      rhl3 = kp*kp+k*k
      rhl4 = 2.*kp*k
c     partial wave projection
      Do 190 ncos=1,ncoss
         x = xcos(ncos)
         q2 = abs(rhl3-rhl4*x)
c     check for rho**2
         If (nifty(16).eq.2) GoTo 100
         If (na.eq.3.and.nifty(7).eq.4) GoTo 130
         If ((q2*b2/8.).ge.150.) GoTo 190
         rho2 = rho2f(b2,alph)
         GoTo 110
c      special calc for ca or other ehere Do rho**2 numercllly
 100     nff = 4
         Call ffact (q2,ff,nff)
         rho2 = ff(3)
c     calc pin com momentun
 110     If (nifty(7).eq.4) GoTo 120
         son = softh(k0,k0)
         sin = softh(k,kp)
         sout = softh(kp,k)
         kap02 = kcm2(son)
         kapp2 = kcm2(sout)
         kap2 = kcm2(sin)
         kap0 = Sqrt(kap02)
         kap = Sqrt(kap2)
         kapp = Sqrt(kapp2)
         facts = gs(kapp)*gs(kap)/gs(kap0)**2
         factp = gp(kapp)*gp(kap)*k0*k0*x/gp(kap0)**2
         revofq(1) = (reb*facts+rec*factp)*rho2*anorm
         imvofq(1) = (imb*facts+imc*factp)*rho2*anorm
         revofq(2) = 0.
         imvofq(2) = 0.
         GoTo 160
c     isobar calc of absp for non he3 nucleus
c     non-a=3 nucleus
c
c     code for tofth(j) for absp amplitudes
c               1,2...pp noflip re,im    3,4...pn noflip re,im
c               5,6...pp flip re,im      7,8...pn flip re,im
 120     rhl6 = rho2*anorm/rhl5
         rhopp = pairpp*rhl6/2.
         rhonn = pairnn*rhl6/2.
         rhopn = pairpn*rhl6
         GoTo 140
c     special 3he calc with isobar absp
 130     rhopp = 0.
         rhonn = 0.
         rhopn = 0.
         kbypp = q2*bp2/8.
         kbynn = q2*bbn2/8.
         kbypn = q2*bp2*bbn2/(bp2+bbn2)/4.
         If (kbypp.lt.150.) rhopp = arhopp*rho2f(bp2,alphp)
         If (kbynn.lt.150.) rhonn = arhonn*rho2f(bbn2,alphn)
         If (kbypn.lt.150.) rhopn = arhopn*exp(-kbypn)*(1.+0.25*alphp*(6
     1   .*bpn2-bpn2*bpn2*q2))
 140     Call tnucth (kp,k,k0,x,tofth,md)
         Do 150 j=1,jmax
            j4 = 4*(j-1)
            revofq(j) = (tofth(1+j4)*(rhopp+rhonn)+tofth(3+j4)*rhopn)
            imvofq(j) = (tofth(2+j4)*(rhopp+rhonn)+tofth(4+j4)*rhopn)
 150     Continue
 160     pl(1) = wt(ncos)
         pl(2) = x*wt(ncos)
         Do 180 j=1,jmax
            revl(1,j) = revl(1,j)+pl(1)*revofq(j)
            imvl(1,j) = imvl(1,j)+pl(1)*imvofq(j)
            revl(2,j) = revl(2,j)+pl(2)*revofq(j)
            imvl(2,j) = imvl(2,j)+pl(2)*imvofq(j)
            nn = nls-1
            If (nn.lt.2) GoTo 180
            Do 170 n=2,nn
               g = x*pl(n)
               pl(n+1) = g-pl(n-1)+g-(g-pl(n-1))/n
               revl(n+1,j) = revl(n+1,j)+pl(n+1)*revofq(j)
               imvl(n+1,j) = imvl(n+1,j)+pl(n+1)*imvofq(j)
 170        Continue
 180     Continue
 190  Continue
      If ((k.eq.kp).and.(k.eq.k0)) Write (6,210) k,imb,reb,imc,rec,(revl
     1(n,1),imvl(n,1),revl(n,2),imvl(n,2),n=1,nls)
      Return
c
 200  format (17h n,sig(mb)p,sigs=,i5,2e15.3)
 210  format (43h true 2n annihilation params for k(=kp=k0):,
     1 50h imb,reb,imc,rec/ then, re, im vabs(l) flip+noflip,/,
     2 f10.3,4e15.5,/,30(4e15.5/))
      End