ci* @(#)xsects.f	1.4      latest revision 6/14/90 17:00:07 */
      subroutine xsects( theta, dsig, sigbrn, dnof, dflip, polar,
     $                   scoul, alph, xgam, sig0, ldmax, noangs,
     $                   nspin, nsex)
********************************************************************************
c     calculate diff xsectns, mod 6/89 by rhl for f amplitude

      implicit real*8 ( a-h, m , o-z )

      character   SCCSID

      complex*16 zi , zcul , zifase
      complex*16 za , zb , zc , zd , ze, zf
      complex*16 zan , zbn , zcn , zdn , zen, zfn
      complex*16 zas(361), zbs(361), zcs(361), zds(361), zes(361),
     &           zfs(361)

      real*8 k0, k02, kone, ktwo, kolmo

      dimension tbr(0:99,12), tbi(0:99,12), tr(0:99,12),ti(0:99,12)
      dimension Pl(0:99), Plp(0:99)
      dimension sigbrn(361) , scoul( 100 )
      dimension sefz(25)
      dimension theta(361), ff(4), dflip(361), dnof(361)
      dimension dsig(361), polar(361)
      dimension nifty(20)

      common /kinemt/   El, Ecm, k0, plab, s , Mp2, Mnuc, Mnuc2
      common /params/   hbarc, pi, Mp, MN, nz, na, nes, nwaves
      common /switch/   nifty
      common /sizes/    achp, acmp, wsp, achn, acmn, wsn
      common /ranges/   rcoul, rcut
      common /inputs/   Tlab , b , ymin1 , ymin2 , xnang,
     $                  kode , lxmax , ngp , nr

      common /Rcomn/    rr , ri , rrb , rib
      common /spins/    nspina , nspinb , nspinc , nspind , nspine
      common /Tcomn/    tr , ti , tbr , tbi


********************************************************************************

c-----calculate angular distribution ---------------------------------------

      zi = ( 0. , 1. )
      k02 = k0 * k0
**************** tvm 8/26/92 reduce output for c13 ****************
      if (na.eq.13) then
      noangs = (60./xnang) + 1
      else
      noangs = (180./xnang) + 1
      endif
********************************************************************
c     write(6,999)xnang,noangs
c999  format('  in xsectns, xnang, noangs=',f5.3,4x, i4)
      cosold = 2.

      do 520 iang = 1,noangs
         theta(iang) = (iang - 1) * xnang
         x = cos(theta(iang) * pi/180.)

c *** rhl chanmge of paez

         if (nifty(6) .ne. 8) then

            call legpol (x, Pl, ldmax)

            gbr = 0.e00
            gbi = 0.e00
            gr =  0.e00
            gi =  0.e00
            if (nifty(6) .eq. 3) then
               call PLPRME (x, Plp, ldmax)
            endif
            sumbr = 0.
            sumbi = 0.
            sumr = 0.
            sumi = 0.

            do 450 ldum = 0, ldmax-1
               sag = ldum + 1
               sumbr = sumbr + (sag * tbr(ldum,1) +
     $                              ldum * tbr(ldum,2)) * Pl(ldum)
               sumbi = sumbi + (sag * tbi(ldum,1) +
     $                              ldum * tbi(ldum,2)) * Pl(ldum)
               sumr = sumr   + (sag * tr(ldum,1)  +
     $                              ldum * tr(ldum,2))  * Pl(ldum)
               sumi = sumi   + (sag * ti(ldum,1)  +
     $                              ldum * ti(ldum,2))  * Pl(ldum)
               if (nifty(6) .eq. 3) then
                  gbr = gbr + (tbr(ldum,1) - tbr(ldum,2)) * Plp(ldum)
                  gbi = gbi + (tbi(ldum,1) - tbi(ldum,2)) * Plp(ldum)
                  gr  = gr  + (tr(ldum,1)  - tr(ldum,2))  * Plp(ldum)
                  gi  = gi  + (ti(ldum,1)  - ti(ldum,2))  * Plp(ldum)
               endif
 450        continue
            rhl = (sumr * hbarc/k0)
            rhl5 = (sumi * hbarc/k0)
            rhl6 = 0.
            rhl7 = (hbarc/k0)**2
            dflip(iang) = rhl7 * (1. - x * x) * (gr**2 + gi**2)
            dflipb = rhl7 * (1. - x * x) * (gbr**2 + gbi**2)
            if (na .eq. 2) dflipb = dflipb * 2./3.
            if (na .eq. 2) dflip(iang) = dflip(iang) * 2./3.
            dnof(iang) = rhl7 * (sumr**2 + sumi**2)
            dnofb = rhl7 * (sumbr**2 + sumbi**2)
            dsig(iang) = dflip(iang) + dnof(iang)
            dsigma = dsig(iang)
            dsigb = dflipb + dnofb

c           *** polarization

            if (dsig(iang) .ne. 0) then
               polar(iang) = 2. * rhl7 * sqrt(1. - x * x) *
     $                        (sumi*gr-sumr*gi)/dsig(iang)
               r = (x * (dnof(iang) - dflip(iang)) - 2. * (1. - x * x) *
     $                  rhl7 * (sumr * gr + sumi * gi))/dsig(iang)
            endif

         else
c           *** use full spin 1/2 x 1/2 amplitude
c           *** assign cross sections for 1/2*1/2

            can1 = k0

            call EXPERA(x,can1,ldmax,za,zb,zc,zd,ze,zf,tr,ti,sefz,scoul)

***************    Temporary write of AMPS to tape 11  *******
*************** T-S mixing 9/19/92 ************************
************* tvm amps 9/19/91 ************************
c       chuck = acos(x) * 180.0 /pi
c       write(11,1966) chuck, za, zb, zc, zd, ze, zf
c1966  format( f6.2, 6(e10.4,e10.4))
************ ******************** **********************

            zan = za
            zbn = zb
            zcn = zc
            zdn = zd
            zen = ze
            zfn = zf
            zas(iang) = za
            zbs(iang) = zb
            zcs(iang) = zc
            zds(iang) = zd
            zes(iang) = ze
            zfs(iang) = zf
            rhl7 = (hbarc)**2

c           *** sig=(a2+b2+c2+d2+e2+f2)2=ds/dw     sig*polar=(a2-b2-c2+d2+e2)/2
c           *** dnof=/a+b/**2     dflip=e2/4=spin-orbit
c           *** pn00, dl0m0, clm00 are various polarization correlatoions

            dsig(iang) = rhl7 * sefz(1)
            dsi0ma = dsig(iang)
            dsigma = dsi0ma
            if (dsigma .eq. 0.) dsigma = 1.0e-10
            polar(iang) =  sefz(3)
            dnof(iang) = rhl7 * sefz(11)
            dflip(iang) = rhl7 * sefz(10)
c ***
c           ***    now do born

            call EXPERA (  x, can1, ldmax, zan, zbn,
     $          zcn, zdn, zen, zfn, tbr, tbi, sefz,scoul)

            dsigb = rhl7 * sefz(1)
         endif

 453     coscm = x
         rhl7 = (hbarc)**2
         tivt = +2. * (k0**2) * (1. - coscm)
         qtf = 0.
         q2f = tivt/rhl7
         if (q2f .ge. 0.) qtf = sqrt(q2f)

c include coulomb calculation with gamma = 0.0
c        if (nifty(10) .eq. 0) go to 510

c *** include coulomb ampl + phase changes

         ff(1) = 1
         if (nifty(10) .ge. 3) go to 500

c ***     reset raddii in case equiv wiped them out

         if (na .eq. 2) go to 480
c        if ((na .le. 16) .AND. (na .gt. 4)) go to 490
c        if (nifty(16) .eq. 0) go to 490

c ****************** ffact now called for any nucleus **********************
c *** for ca or non-gaussian nucleus

         call ffact (+tivt, ff, 1)

         go to 500
 480     qr = qtf * rcoul
         if (qr .eq. 0.) go to 500
         ff(1) = (3./qr/qr) * (sin(qr)/qr - cos(qr))
         go to 500
 490     ff(1) = (1. - alph * (qtf * acmp/hbarc)**2) *
     $            exp(-.25 * (qtf * achp/hbarc)**2)

c *** the value of theta(iang) from e-10 to e-06 necessary for vax

 500     if (theta(iang) .eq. 0) theta(iang) = 1.e-06
         phasec = sin(theta(iang) * pi/360.)
         fcoul = (-xgam * ff(1)/2./k0/phasec/phasec) * hbarc
         dsigcl = fcoul * fcoul * 10.
         phasec = -xgam * LOG(phasec**2)
c        tm 9/19/90 + 2. * sig0 removed in next two lines
         rhl6 = phasec + 2.0 * sig0
         phasec = phasec  + 2.0 * sig0
         zifase = zi * phasec
         zcul = fcoul * (cos(phasec) + zi * sin(phasec))
         if (nifty(6) .ne. 8.) go to 501

c *** spin 1/2 x 1/2
c *** include coulomb for proton-nucleus by calling  the
c *** next subroutine which adds to za and zb the coulomb
c *** amplitudes. zc ad ac are unnafected.see stapp yps. metr.

********* Halloween only change, tvm 10/30/92 ******************
        index = 1
        nth = 1
c       call yennie (nth,k0,xgam,scoul,coefr,coefi,index,ldmax,x,zcul)
*****************************************************************
         call EXCULA (x,can1,ldmax,za,zb,zc,zd,ze,zf,zcul,sefz)

            zas(iang) = za
            zbs(iang) = zb
            zcs(iang) = zc
            zds(iang) = zd
            zes(iang) = ze
            zfs(iang) = zf
         dsig(iang) = rhl7 * sefz(1)
         dsigma = dsig(iang)
***************    Temporary write of amps to tape 11  *******
************* tm amps 9/19/91 ************************
******************** 10/23/92 this is the one you want tvm write(11 *
        chuck = acos(x) * 180.0 /pi
c       write(11,1966) chuck, za, zb, zc, zd, ze, zf
 1966   format(f6.2, 6(e10.4,e10.4))
************ ******************** **********************
         go to 510
c *** regular cases of either 0 x 0 or 0 x 1/2

 501     rhl = (sumr * hbarc/k0) + fcoul * cos(phasec)
         rhl5 = (sumi * hbarc/k0) + fcoul * sin(phasec)
         dsigma = rhl * rhl + rhl5 * rhl5 + dflip(iang)
         dsi0ma = dsigma
         rhl = (sumr * hbarc/k0)
         rhl5 = (sumi * hbarc/k0)
         sumbr = sumbr * hbarc/k0 + fcoul * cos(rhl6)
         sumbi = sumbi * hbarc/k0 + fcoul * sin(rhl6)
         rhl6 = fcoul/sqrt(dsig(iang))
         dsigb = sumbr * sumbr + sumbi * sumbi + dflipb
         dnof(iang) = rhl * rhl + rhl5 * rhl5
         dnofb = sumbr**2 + sumbi**2

c *** calc of t.f. to lab from cm, see goldberg+watson
c ***

c move total cross section calculation out of IANG do loop
 510     continue
         ektwo = Mnuc + plab * plab * Mnuc * (1. - coscm)/s
         ktwo = sqrt(ektwo**2 - Mnuc2)
         ekone = el + Mnuc - ektwo
         kone = sqrt(abs(ekone**2 - Mp2))
         coslab = sqrt(abs(kone**2 - k02 * (1. - coscm**2)))/kone
         if (coslab .gt. cosold) coslab = -coslab
         cosold = abs(coslab)
         dot = (Mnuc2 + Mp2 + 2. * ekone * ektwo - s)/2.
         cmtol = kone * kone * s/
     $                  (plab * Mnuc *(ektwo * kone - dot * ekone/kone))

c *** convert cross sections from fm**2 to mb

         if (nifty(6) .eq. 3) dsi0ma = dsigma
         dsi0ma = dsi0ma * 10.
 513     dnof(iang) = dnof(iang) * 10.
         thlab = acos(coslab) * 180./pi
         dsigma = dsigma * 10.
         dsigb = dsigb * 10.
         dsig(iang) = dsig(iang) * 10.
         sigbrn(iang) = dsigb
         dnofb = dnofb * 10.
         dflip(iang) = dflip(iang) * 10.
         dflipb = dflipb * 10.
         dsigl = dsigma * cmtol
         dsigbl = dsigb * cmtol

c *** print out of cross sections vs theta

         if (nifty(6) .ne. 8) then
            write(6,1100) theta(iang), coscm, dsi0ma, dsigb,
     $                    q2f, thlab, coslab, dsigl, dsigbl,
     $                    rhl, rhl5, dsig(iang), dsigcl
            if (nifty(20) .gt. 1) then
c              write(10,4000) theta(iang), dsi0ma, polar(iang), r
            endif
         endif
         if (nifty(6) .eq. 3) then
            write(6,1120) dnof(iang), dnofb, polar(iang),
     $                    dflip(iang), dflipb, r
         endif
 1811    if (nifty(6) .ne. 8) go to 1986

c *** spin 1/2 * 1/2 print out

c	mod 6/90 for sefz(15) for 2' printout
*********** tvm 3/8/92 temp modification for 800 MeV He3 **********
c        queeg = 180./pi*dasin(955.0/1463.*sin(pi*theta(iang)/180.))
c        write(6,6510) queeg, qtf, dsig(iang), polar(iang),
c    $                 sefz(6),sefz(7),sefz(8),sefz(9),sefz(2),sefz(15)
********************  temp write above *********************************
         write(6,6510) theta(iang), qtf, dsig(iang), polar(iang),
     $                 sefz(6),sefz(7),sefz(8),sefz(9),sefz(2),sefz(15)
c    **************    Ratio stuff *****************
       if(nz.eq.1.and.nifty(1).eq.5) then
       write(7,7844) theta(iang), dsig(iang)
       endif
       if(nz.eq.2.and.nifty(1).eq.5) then
       write(8,7844)  theta(iang), dsig(iang)
       endif
       if(nz.eq.1.and.nifty(1).eq.9) then
       write(88,7844)  theta(iang), dsig(iang)
       endif
       if(nz.eq.2.and.nifty(1).eq.9) then
       write(77,7844)  theta(iang), dsig(iang)
       endif
 7844  format(2e14.4)
***** tvm 3/2/92  temp write  *****************************
c      solo = - sefz(7) 
c      tvm = dsig(iang) * 3.14159265 * 1000000.0 / (k0 * k0)
c      write(11,6511) queeg, tvm, solo
c6511    format(1h , f6.2, 2e12.4)
************************ end temp modifications *************
 6510    format(1h ,f6.2,f8.4,8e12.4)
c      write(7,7311) theta(iang), sefz(20), sefz(21), sefz(22), sefz(23)
c7311  format(1h , f6.2, 4e12.4)
         if (nifty(20) .gt. 1) then
            write(10,4010) theta(iang), dsig(iang), sefz(3),sefz(2),
     $                     sefz(15),sefz(6),sefz(7),sefz(8)
c	plot variables on tape10
c      	th,sig,#2-Aoono,#11-Aoonn,#2'-Aooon,#3-Dnn,#4-Dsl,#5-Dss
         endif
c        write(6,1120) za, zb, zc, zd, ze, zf
 1986    dsig(iang) = dsigma
         if (nifty(6) .eq. 3 ) go to 520

c *** set up for plot, polar for 1/2x1/2
         dnof(iang) = sefz(2)
         dflip(iang) = sefz(15)
         dsig(iang) = dsigma
 520  continue

c *** Calculate total cross-sections or store t's for Cex ***

c     call TOTSIG( k0, scoul, sigl, xgam, ldmax, nspin, nsex )

         if (nifty(10) .eq. 0) then
c total cross section calculation moved out of IANG do loop
            xrhl=1.0
            call EXPERA( xrhl, can1, ldmax, za, zb,
     $               zc, zd, ze, zf, tr, ti, sefz)
         endif

      rhl = 10. * (hbarc)**2
      sefz(16) = sefz(16) * rhl
      sefz(17) = sefz(17) * rhl
      sefz(18) = sefz(18) * rhl
      if (nifty(6) .eq. 8) then
         write(6,1921) sefz(16), sefz(17), sefz(18), Tlab, xgam
      endif

      if (nifty(6) .eq. 8) then
      write(6,6540)
 6540 format(1h ,'th-c.m.  ',11x,'za',16x,'zb',16x,'zc',16x,
     $ 'zd',16x,'ze',16x,'zf')
      do 550 iang = 1,noangs
         theta(iang) = (iang - 1) * xnang
         coscm = cos(theta(iang) * pi/180.)
         qtf = 0.
         q2f = 2.*k0*k0*(1.-coscm)/hbarc/hbarc
         if (q2f .ge. 0.) qtf = sqrt(q2f)
*********** tvm 9/26/92 Remove printing of complex a-f amps *******
c        write(6,6550)theta(iang),zas(iang),zbs(iang),zcs(iang),
c    $ zds(iang),zes(iang),zfs(iang)
       toy = 197.32 / sqrt( dsig(iang) /10)
       tva = zas(iang) * toy
       tvb = zbs(iang) * toy 
       tvc = zcs(iang) * toy
       tvd = zds(iang) * toy
       tve = zes(iang) * toy 
       tvf = zfs(iang) * toy
       tma = -zi * zas(iang) * toy
       tmb = -zi * zbs(iang) * toy 
       tmc = -zi * zcs(iang) * toy
       tmd = -zi * zds(iang) * toy
       tme = -zi * zes(iang) * toy
       tmf = -zi * zfs(iang) * toy
       write(11,6668) theta(iang), tva, tma, tvb, tmb, tvc, tmc,
     $                tvd, tmd, tve, tme, tvf, tmf
 6550    format(1h ,f6.2,1x,6('(',e8.2,',',e8.2,')'))
 6555    format(1h ,f6.2, 12e8.2)
 6668    format(f8.2, 12e9.2)
 550     continue
      endif

      RETURN

c *** FORMATS ***

 1100 format(1h ,f7.1,f8.3,2e11.3,2f8.1,f10.3,2e11.3,2f9.3,e12.4
     1,e10.2)
 1120 format (1h ,18x,10e11.3)
 1921 format(1h ,' sigot= ',e13.5,' sigit= ',e13.5,
     $           ' sig2t= ',e13.5,' e= ',f8.2,' xgam= ',e13.5)
 4000 format(' ', f7.2, 3e11.3)
 4010 format(1h ,f7.2,2e11.3)
       end
