subroutine tnucth (jj1,jj2,kp,k,k0,cthnuc,tofth,mn)
c *** see r h 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)) go to 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.) go to 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) go to 70
      if (nifty(5).eq.9) go to 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) go to 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) go to 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) go to 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