c/* @(#)thep.f	1.2      latest revision 6/26/89 09:07:59 */
      subroutine thep(xx,k0,neles,za,zb,zc,zd,ze,zf,tr,ti,scoul)
************************************************************************

c     forms the p-He3 amplitudes using the formulas is Stapp
c     rhl mod 6/89 for possible z ampliutde; YET f set =0 here
c     RHL mod 4/83 to mult nspin=1,2 by 2i to agree with Stapp
c     Ypsilantis and Metropolis Phys. Rev. 105 1957

      implicit real*8 (a-h, k, o-y)
      implicit complex*16 (z)

      dimension y(100), y1(100), y2(100), a(5,5)
      dimension ti(100,12), tr(100,12)
      dimension scoul(100)

c    optical potential values passed via common/vopt/

      common /vopt/ v10(2), v1m1(2), v01(2), v11(2), v00(2), vss(2),
     &              v1s(2),vsum1, vsum2

c >>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

      zum1 = (0.0, 0.0)
      zum2 = (0.0, 0.0)
      zum3 = (0.0, 0.0)
      zum4 = (0.0, 0.0)
      zum5 = (0.0, 0.0)
      zum6 = (0.0, 0.0)
      zum7 = (0.0, 0.0)
      zum8 = (0.0, 0.0)
      zum9 = (0.0, 0.0)

      call legpol(xx,y,neles)
      call plprme(xx,y1,neles)
      call pldblp(xx,y2,neles)

      zi = (0.0, 1.0)
      seno=+sqrt(1.-xx*xx)
      seno2=seno*seno
      do 10 l=1,neles
      call costap(a,l)
c
c     define amps as in stapp paper with
c          zts=z6=singlet  z1=t11     z2=t00     z3=t01
c          z4=t10         z5=t1-1    z7=t01/sinth
****             z9 -> triplet-singlet mixing
c      t1,2,3,4,5 refer to stapp convwention whic ne nspin(rhl)
c          nstapp-nspin(rhl)
c        2-2   1-3  5-4  3-5  4-6   s-1
      zts=(tr(l,1)+zi*ti(l,1))*2.*zi
      zt1=(tr(l,3)+zi*ti(l,3))*2. * zi
      zt2=(tr(l,2)+zi*ti(l,2))*2.*zi
      zt4=(tr(l,6)+zi*ti(l,6)) * 2. * zi
      zt9 = 0.
      if(l.gt.1) zt9=(tr(l,7)+zi*ti(l,7))*2.*zi
      if(l.lt.3)go to 20
c     l=ldum ge 3, no special cases
      Lm2=l-2
      zt3=(tr(Lm2,5)+zi*ti(Lm2,5)) * 2. * zi
      zt5=(tr(Lm2,4)+zi*ti(Lm2,4)) * 2. * zi
      go to 30
 20   zt5=0.
      if(l.eq.2)go to 25
c     l=1 (lactual=0)
      zt2=0.
      zt3=0.
      go to 30
c     l=2 (lactual=1), special use of r11**0 now done here
 25   zt3=(tr(1,2)+zi*ti(1,2) )*2.*zi
 30   continue
c     ** zz factor is correction for inclusion of e^2i(sig) in *all* channels
      zz = zexp(-2.0 * zi * scoul(l))
      z1=a(1,1)*zt1+a(1,2)*zt2+a(1,3)*zt3+a(1,4)*zt4+a(1,5)*zt5
      z2=a(2,1)*zt1+a(2,2)*zt2+a(2,3)*zt3+a(2,4)*zt4+a(2,5)*zt5
      z3=a(3,1)*zt1*zz+a(3,2)*zt2*zz+a(3,3)*zt3*zz+a(3,4)*zt4+a(3,5)*zt5
      z7=a(3,1)*zt1*zz+a(3,2)*zt2*zz+a(3,3)*zt3*zz+a(3,4)*zt4+a(3,5)*zt5
      z4=a(4,1)*zt1*zz+a(4,2)*zt2*zz+a(4,3)*zt3*zz+a(4,4)*zt4+a(4,5)*zt5
      z8=a(4,1)*zt1*zz+a(4,2)*zt2*zz+a(4,3)*zt3*zz+a(4,4)*zt4+a(4,5)*zt5
      z5=a(5,1)*zt1*zz+a(5,2)*zt2*zz+a(5,3)*zt3*zz+a(5,4)*zt4+a(5,5)*zt5
      al=l-1
      z6=0.5*(2.*al+1.)*zts
      z9 = 0.0
      if (l.gt.1) z9=-zt9 * sqrt(2.)*(2.*al+1.)/sqrt(al*(al+1.))/4.
c     multiply by appropiate leg pol
      z9=-zi*z9*seno*y1(l)
      z6=-zi*z6*y(l)
      z1=-zi*z1*y(l)
      z2=-zi*z2*y(l)
      z3=-zi*z3*seno*y1(l)
      z4=-zi*z4*seno*y1(l)
      z5=-zi*z5*seno2*y2(l)
c     leave out sin(theta) in 7 and 8 to avoid /0 later
       z7=-zi*z7*y1(l)
       z8=-zi*z8*y1(l)
      zum1=zum1+z1
      zum2=zum2+z2
      zum3=zum3+z3
      zum4=zum4+z4
      zum5=zum5+z5
      zum6=zum6+z6
      zum7=zum7+z7
      zum8=zum8+z8
      zum9=zum9+z9
 10   continue
      zss=zum6/k0
      z11=zum1/k0
      z00=zum2/k0
      z01=zum3/k0
      z10=zum4/k0
      z1m1=zum5/k0
      z1s=zum9/k0
c special for d
      z01es=zum7/k0
      z10es=zum8/k0
****      tvm 10/5/92   write amps to tape7 as a test     *****
       pi = 3.14159
       th = acos(xx) * 180./pi
      pss = zss
      p00 = z00
      p11 = z11
      p10 = z10
      p1s = z1s
      qss = -zi * zss
      q00 = -zi * z00
      q11 = -zi * z11
      q10 = -zi * z10
      q1s = -zi * z1s
c     write(7,333) th, pss, qss, p00,q00,p11,q11,p10,q10, p1s, q1s
 333   format(f5.1,10e12.4)

c     form a,b,c,d,e f of bistricky,lehar,winternitz
      za=0.5*(z11+z00-z1m1)
      zb=0.5*(z11+zss+z1m1)
      zc=0.5*(z11-zss+z1m1)
c      rhl change 4/83 include - sign
      zd=-(z10es+z01es)/sqrt(2.)
c    ***** tm amps 06/9/91 insert - sign in e ******
c     ze = -zi * (z10-z01)/sqrt(2.)
      ze= zi * (z10-z01)/sqrt(2.)
      zf = -zi * sqrt(2.) * z1s

      call voptth(k0, k0, k0, cthnuc)
      return
c**  this program valid on ftn4 and ftn5 **
      end