Subroutine tnucth (jj1,jj2,kp,k,k0,cthnuc,tofth,mn)
c *** see RH Landau, Program lpott, 1981
c     tnucth calculates the theta dependent ampl in pi-nucleus frame
c     with angle dependent po
      Implicit Real*8 (a-h,i,k,m,o-z)
      Dimension nifty(20), bnuc(20), bnucf(20), bn(16), bnf(16), retij
     1(14), imtij(14), tofth(8)
      Dimension km(4), qim(4), em(4), qm(2), pm(4), kapm(4), alph(4),
     1 beta(4)
      Common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,xmn,nz,
     1nes,nwaves,nifty,na
c     tofth(j) code, t = ...,for j=...
c     1    retp-spinindep              2    imtp-spin indep
c     3    retn- spin indep            4    imtn-spin indep
c     5    retp-spin dependt           6 imtp-spin dep
c     7    retn-soin dept              8    imtn spon-depend
c     statement functions now follow
      epi(p) = Sqrt(mpi2+p*p)
      en(p) = Sqrt(mn2+p*p)
c     the pi-nucleon c.o.m.. momemtum**2 for given s
      kcm2(s) = (s-(mpi+mn)**2)*(s-(mpi-mn)**2)/4./s
c     internal nucleon energy, new optimal definition
      eint(p,pp) = Sqrt(mn2+(rhl1**2)*p*p*0.25+(rhl2**2)*pp*pp*0.25-0.5*
     1rhl1*rhl2*p*pp*cthnuc)
c     softh is the new,optimal,theta-dependent s
      softh(p,pp) = mpi2+mn2+rhl1*p*p-rhl2*p*pp*cthnuc+2.*epi(p)*eint(p,
     1pp)
      nit = 1
      amass = na*mn/xmn
      ma = na*xmn
      rhl1 = (amass+1.)/amass
      rhl2 = (amass-1.)/amass
      mn2 = mn*mn
      mpi2 = mpi*mpi
      m12 = (ma-mn)**2
      rhl3 = mpi2-mn2
      ma2 = ma*ma
      sin = softh(k,kp)
      sout = softh(kp,k)
      enin = eint(k,kp)
      enout = eint(kp,k)
      If ((nifty(5).gt.5).and.(nifty(5).ne.9)) GoTo 20
c     calculate on-shell kappa with 2 body energy
      son = softh(k0,k0)
      kapon2 = kcm2(son)
      kapon = Sqrt(kapon2)
      kap2 = kcm2(sin)
      kap = Sqrt(kap2)
      kapp2 = kcm2(sout)
      kapp = Sqrt(kapp2)
      If (kapp.le.0.) kapp = 0.0001
      If (kap.le.0.) kap = 0.0001
      cthn = ((epi(kap)*epi(kapp)-epi(k)*epi(kp))/kapp/kap)+k*kp*cthnuc/
     1kap/kapp
      If (nifty(5).eq.2) cthn = cthnuc
      If (abs(cthn).le.2.) GoTo 10
      Write (6,100) cthn
 10   Continue
      gamth = Sqrt(epi(kap)*epi(kapp)*en(kap)*en(kapp)/(epi(k)*epi(kp)*
     1enin*enout))
      costh = cthn
      factor = gamth
c     choose defintion of energy
c     nifty(5)=0   % use softh(on shell k)
c              1 use sin of theata (incoming momentum)
c              9 use softh(on shell k)plus aay angles
c              2 softh (on shell), no angle t.f.
      kape = kapon
      If (nifty(5).eq.1) kape = kap
      ekape = epi(kape)+en(kape)
      If (nifty(5).lt.5) GoTo 70
      If (nifty(5).eq.9) GoTo 40
c     3 body def of energy
c     nifty(5)= 5: e3b with nonfolded or w0-folded t + old,optimal angle
c             =7: same as above yet with aay magic angle t.f.(standard)
c     n.b. n(5)=7 is the standard case, whereas 6 is for trial
c     Dont use 6 If have w0-folded t s
 20   Continue
      spinuc = epi(k0)+Sqrt(ma2+k0*k0)
      q2 = kp*kp+k*k-2.*kp*k*cthnuc
      qk = k*kp*cthnuc-k*k
      be = 5.*nifty(2)
      pkap2 = rhl2*rhl2*(k*k+q2/4.+174.*174.+qk)
      ekape = spinuc-be-Sqrt(pkap2+m12)
      ekape = ekape*ekape-pkap2
c     e3b-folded t s being used
c       non reltv def for folding
      If (nifty(5).ne.6) GoTo 30
      twomu = 2.*(mpi+mn)*(ma-mn)/(mpi+ma)
      ekape = spinuc-ma+mn+5-be-(k*k+(q2-k0*k0)/4.+qk)*rhl2*rhl2/twomu
 30   rhl = 1.
      If (ekape.lt.0.) rhl = -1.
      ekape = rhl*Sqrt(rhl*ekape)
      If (nifty(5).eq.5) GoTo 70
c     nifty(5) = 6,7,9 use magic vector prescription(aay)
c     1st two components x,y of kin, 3rd 4th,x y of kout
 40   rhl = Sqrt(1.-cthnuc**2)
      ki0 = epi(k)+eint(k,kp)
      kf0 = epi(kp)+eint(kp,k)
      km(1) = k
      km(2) = 0.
      km(3) = kp*cthnuc
      km(4) = kp*rhl
      Do 50 j=1,2
         j2 = j+2
         qm(j) = km(j2)-km(j)
         pm(j) = -km(j)/amass+rhl2*qm(j)/2.
         pm(j2) = pm(j)-qm(j)
         kapm(j) = km(j)+pm(j)
         kapm(j2) = kapm(j)
         alph(j) = rhl3/sin
         alph(j2) = rhl3/sout
         qim(j) = 0.5*(km(j)-pm(j)-alph(j)*kapm(j))
         qim(j2) = 0.5*(km(j2)-pm(j2)-alph(j2)*kapm(j2))
         beta(j) = qim(j)*kapm(j)/ki0/(ki0+Sqrt(sin))
         beta(j2) = qim(j2)*kapm(j2)/kf0/(kf0+Sqrt(sout))
         If (j.ne.2) GoTo 50
         beta(1) = beta(1)+beta(2)
         beta(2) = beta(1)
         beta(3) = beta(3)+beta(4)
         beta(4) = beta(3)
 50   Continue
      Do 60 j=1,4
         em(j) = qim(j)-beta(j)*kapm(j)
 60   Continue
      kap = Sqrt(em(1)**2+em(2)**2)
      kapp = Sqrt(em(3)**2+em(4)**2)
      costh = (em(1)*em(3)+em(2)*em(4))/kap/kapp
      factor = Sqrt(epi(kap)*epi(kapp)*en(kap)*en(kapp)/(epi(k)*epi(kp)*
     1enin*enout))
c     run with shifted pi n subenergy( for cex)
 70   eshift = ekape
      If (nifty(1).eq.1) eshift = eshift+nifty(11)
      If (nifty(1).eq.0.or.nifty(1).eq.8) eshift = eshift+nifty(12)
      If (nifty(1).eq.4.or.nifty(1).eq.7) eshift = eshift+nifty(11)
      If (mn.lt.1000.) Call tnoff (jj1,jj2,kapp,kap,eshift,kp,k,cthnuc)
c       deactivate call for now
c     ----  If(mn.gt.1000.)Call tDoff(kapp,kap,eshift,kp,k,cthnuc)
      x = costh
      pl2 = 1.5*x*x-.5
      pl3 = 2.5*x**3-1.5*x
c     indepentent amplitudes
      Do 80 j=1,4
         tofth(j) = factor*(bn(j)+x*bn(j+4)+pl2*bn(j+8)+pl3*bn(j+12))
 80   Continue
      If (nifty(6).ne.3) Return
c     dependent amplitudes
      rhl = factor*Sqrt(abs(1.-x*x))
      If (abs(x).gt.1.) rhl = -rhl
      pl2 = 3.*x
      pl3 = 0.5*(15.*x*x-3.)
      Do 90 j=1,4
         tofth(j+4) = rhl*(bnf(j+4)+bnf(j+8)*pl2+bnf(j+12)*pl3)
 90   Continue
      Return
c
 100  format (1x,' cos thea pi-n gt 2 = ',e15.3)
      End