c$debug
c***********************************************************************
      SUBROUTINE XSECT
c
c     Subroutine for calculation of cross sections and spin observables
c     from the S matrix elements.
c     
c***********************************************************************
c
C     KMAX  = Number of integration points
C     NCHAN = Number of sub channels in INTEG2
c     ispc  = Labeled common sizes
c     nlev  = Number of levels
c
      parameter (kmax = 300, nchan = 64, ispc0 = 16384, ispc1 = 3*16384
     1, ispc4 = 8192, nlkmax = ispc1/2/kmax, nlev =12)
      parameter (n_output = 7, nmax = 200, nsig = 4, npol = 10)
      logical isym, i_open20, i_out20, i_20flag
c
      COMMON ALPHA(15),HBARC,AMU,CHGSQ,DRT,RMAX
     1,ICON(20),LMAX,KT,MAXD,IOUT,JOUT,NLTR,IDAT(6)
     2,KLAM(nlev),LTRT(nlev),JTRT(nlev),NADD(nlev)
      COMMON EN(nlev),FM(nlev),Z(nlev),FMA(nlev),ZA(nlev),RY(nlev)
     1,FS(nlev),ECM(nlev),FMU(nlev),RC(nlev),FK2(nlev),FK(nlev)
     2,ETA(nlev),AFACT(nlev),PFACT(nlev),DRF(nlev),FN(nlev)
     3,PNLOC(nlev),DUM(nlev),TSIG(nlev)
     4,IS(nlev),NS(nlev),NLI(nlev),IN(nlev),INN(nlev,nlev)
     5,LAM(nlkmax),ISP(nlkmax),JSP(nlkmax),ILK(nlkmax)
     6,LKTAB(2*nlkmax),ICUP(2*nlkmax),JCUP(2*nlkmax)
c
      COMMON /label0/space0(ispc0)
      COMMON /label1/space1(ispc1)
      COMMON /label4/space4(ispc4)
      COMMON /LLL   /NLVAL,NBLKS,NLMIN,NLDEL,MAXTAB(nlev)
     1,INDTAB(nlev+2)
      common /xsectc/thetan,theta1,dtheta,angdat(nmax)
      common/ch3files/input, ioutp
c
      DIMENSION  D(ispc0), Csig(nmax), Polz3(npol)
     1    , Sigma(nmax,n_output+1), PLM(1000), SR(1000)
c
      EQUIVALENCE
     1  (Space0(    1),D(1)),    (Space1(1), sr(1))
     2, (Space1(ispc0+1), PLM(1))
     3, (Space4(1), Sigma(1,1)), (Sigma(1,8), Csig(1))
c
c	Sigma(i , 1) stores  SIG  Cross section
c	Sigma(i , 2) stores  POL  Polarization
c	Sigma(i , 3) stores  ASY  Asymmetry
c	Sigma(i , 4) stores  AYY  
c	Sigma(i , 5) stores  A22  Tensor asymmetry 20
c	Sigma(i , 6) stores  A21  Tensor asymmetry 21
c	Sigma(i , 7) stores  A20  Tensor asymmetry 22
c	Sigma(i , 8) stores  Csig SIG/Rutherford 
c
      DATA PI, i_20flag /3.141592654, .true./
C
      i=0
      mks=2.0*maxtab(1)
      do 30 ii=1,nltr
      maxtab(ii)=maxtab(ii)*mks
      i=i+maxtab(ii)
   30 continue
      if(i .gt. ispc0) then
        write(ioutp,'(2(a,i5))')'0Amplitude storage too small, need',i
     1,' ,have',ispc0
        return
      endif 
c
      if(fm(1).eq.fma(1).and.z(1).eq.za(1).and.is(1).eq.jtrt(1)) then
        isym = .true.
      ELSE
        isym = .false.
      ENDIF
c
      JMAX=LMAX
      LMAX2=LMAX+LMAX
      NTHETA=THETAN
      NTHETA=MIN0(NTHETA,NMAX)
      IF(NTHETA.EQ.0) GO TO 3001
      MAX1=JTRT(1)+1
      MAXM=0
      DO 60 I=1,NLTR
      MAXM=MAX0(MAXM,(JTRT(1)+JTRT(I)+IS(1)+IS(I))/2+1)
   60 CONTINUE
c
      tint=THETA1
      do 80 nth=1,Ntheta
      if(icon(1).eq.3.or.icon(1).eq.7) then
        if    (tint.gt. 1.0) then
          tint = 1.0
        elseif(tint.lt.-1.0) then
          tint =-1.0
        endif
        theta = acos(tint)*180.0/pi
      elseif(icon(1).eq.4) then
        theta=angdat(nth)
      elseif(icon(1).eq.5) then
        theta=acos(angdat(nth))*180./pi
      else
        theta = tint
      endif
      angdat(nth) = theta
      tint = tint + dtheta
   80 continue
C
C     CHANNEL LOOP
C
      DO 3000 I=1,NLTR
      MAXI=JTRT(I)+1
      JR=NS(1)
      JS=NS(I)
      maxii=maxtab(i)
        if(i.eq.1 .and. icon(6).eq.0
     1   .and. iand(icon(17),1).eq.0) go to 1000
C
C     ANGLE LOOP
C
      DO 980 NTH=1,NTHETA
      theta = angdat(nth)
      CALL LGNDR(PLM,MAXM,LMAX,THETA)
c
      INC=INDTAB(I)+1
      DO 110 J=1,maxii
      SR(J)=0.0
  110 CONTINUE
C
      DO 560 NB=1,NBLKS
      JA1=1
      JA2=LMAX
C     READ IN AMPLITUDE BLOCKS
      IF(NBLKS .ne. 1) then
        DO 190 II=1,NLTR
        READ (2)I1,I2,JA1,JA2,(D(IT),IT=I1,I2)
  190   CONTINUE
      ENDIF
c
      DO 500 L=JA1,JA2
      LZ=L-1
      IT=1
      MM=-JTRT(1)
c	initial target
      DO 490 MX=1,MAX1
      MMS=-JTRT(I)
c	final target
      DO 480 MY=1,MAXI
      MS=-IS(1)
c	initial projectile
      DO 470 M =1,JR
      MSP=-IS(I)
c	final projectile
      DO 460 MP=1,JS
      ML=MS+MM-MSP-MMS
      MLT=IABS(ML/2)
      INDX=MLT*LMAX+LZ+1
      SR(IT  )=SR(IT  )+D(INC  )*PLM(INDX)
      SR(IT+1)=SR(IT+1)+D(INC+1)*PLM(INDX)
      INC=INC+2
      IT=IT+2
      MSP=MSP+2
  460 CONTINUE
      MS=MS+2
  470 CONTINUE
      MMS=MMS+2
  480 CONTINUE
      MM=MM+2
  490 CONTINUE
  500 CONTINUE
  560 continue
      rewind 2
c
      IT = 1
      IP = 1
      DO 790 MX=1,MAX1
      DO 780 MY=1,MAXI
C
      IF(I.EQ.1) THEN
      IF(ETA(1).EQ.0.0) THEN
        FFACT=1.0
        gfact=0.0
      ELSE
C
C     CALCULATE COULOMB AMPLITUDE FOR CHANNEL 1
C
      IF(abs(theta).lt.0.0625
     1      .or.(isym.and.abs(theta-180.).lt.0.0625)) then
        FFACT=1.0
        GFACT=0.0
      ELSE
        SS=SIN(THETA/114.59156)**2
        ARG=ETA(1)* LOG(SS)
        FFACT=ETA(1)/(2.0*FK(1)*SS)
        FR=COS(ARG)*FFACT
        FI=SIN(ARG)*FFACT
        if(isym) then
        CC=COS(THETA/114.59156)**2
        ARG=ETA(1)*LOG(CC)
        GFACT=ETA(1)/(2.0*FK(1)*CC)
        GR=COS(ARG)*GFACT
        GI=SIN(ARG)*GFACT
      endif
      ENDIF
      DO 620 M =1,JR
      DO 610 MP=1,JS
C
C     ADD COULOMB TERM TO CHANNEL 1 AMPLITUDE
C
      if(isym) then
c       identical particles
        m1=-is(1)+2*(m -1)
        m2=-is(1)+2*(mx-1)
        m3=-is(1)+2*(mp-1)
        m4=-is(1)+2*(my-1)
        if(m1+m2 .eq. m3+m4) then
        is0=max (0,abs(m1+m2))
        is1=2*is(1)
        do 605 iss=is0,is1,2
        temp = vcc(is(1),is(1),iss,m1,m2)*vcc(is(1),is(1),iss,m3,m4)
        t1=(fr + phasef(iss/2)*gr)*temp
        t2=(fi + phasef(iss/2)*gi)*temp
	if(nth .eq. 19) then
	endif
        sr(IT  )=sr(IT  )-t1
        sr(IT+1)=sr(IT+1)+t2
  605   continue
        endif
      else
c       non-identical particles
        IF(M .EQ.MP .AND. MX.EQ.MY) THEN
          SR(IT  )=SR(IT  )-FR
          SR(IT+1)=SR(IT+1)+FI
        ENDIF
      endif
      IT=IT+2
  610 CONTINUE
  620 CONTINUE
      ENDIF
C
      if(isym) then
        Csig(nth)=ffact**2 + gfact**2
     1    +(2.0*phasef(is(1))/float(jr))*(fr*gr+fi*gi)
      else
        Csig(nth)=FFACT**2
      endif
      ENDIF
c
      IP=IP+2*JR*JS
  780 CONTINUE
  790 CONTINUE
c     --------------------------------------------------------
c       Calculate cross section and spin observables and
c       output to disk file 20
c
        i_out20 = .false.
      if(iand(icon(17),2**(i-1)) .ne. 0) then
        i_out20 = .true.
        if(i_20flag) i_open20 = .true.
        i_20flag = .false.
      endif
c
      call polfct(max1,maxi,jr,js,theta,Polz3,sr(1)
     1   ,i_open20,i_out20,nth,ntheta,ALPHA,IDAT)
c     -------------------------------------------------------------------
cc
      do 800 ii=1,n_output
      Sigma(nth,ii) = Polz3(ii)
  800 continue
  980 CONTINUE
c
 1000 CONTINUE
C
      if(icon(19).ne.0) go to 3000
      if(i.eq.1 .and. icon(6).eq.0) go to 2001
      WRITE(ioutp,9901)'Sigma    ',I,LTRT(I)*JTRT(I),ALPHA
     1	,(IDAT(M),M=1,3)
      IF(I.EQ.1 .and. eta(1).ne.0.0) THEN
        WRITE(ioutp,9027)
      ELSE
        WRITE(ioutp,9028)I
      ENDIF
C
      tint=THETA1
      DO 2000 NTH=1,NTHETA
      theta = angdat(nth)
c
      IF(I.EQ.1) THEN
      IF(ETA(1).NE.0.0 .AND. (abs(theta).lt.0.0625
     1      .or. (isym.and.abs(theta-180.).lt.0.0625))) then

        RATIO=1.0
	do 1120 ii=1,n_output
	Sigma(nth,ii) = 0.0
 1120	continue
        ELSE
        RATIO=Sigma(nth, 1)/Csig(nth)
      ENDIF
        Csig(nth)=RATIO
      ENDIF
 1900 CONTINUE
C
C     PRINT CROSS SECTIONS
C
      CS=COS(THETA*pi/180.)
      IF(I.EQ.1 .and. eta(1).ne.0.0) THEN
        WRITE(ioutp,9030)THETA,Csig(nth)
     1  ,(Sigma(nth,ii),ii=1,n_output),Theta,cs
      ELSE
        WRITE(ioutp,9031)THETA
     1  ,(Sigma(nth,ii),ii=1,n_output),Theta,cs
      ENDIF
      tint = tint + DTHETA
 2000 CONTINUE
 2001 continue
      IF(I.EQ.1) THEN
        WRITE(ioutp,9905)'0Reac_Sig',TSIG(1)
      ELSE
        WRITE(ioutp,9905)'0Totl_Sig',TSIG(I)
      ENDIF
      if((icon(6).eq.0.and.i      .eq.1).or.
     1   (icon(1).eq.4.or. icon(1).eq.5))   go to 3000
      IF(ICON(9).NE.0) THEN
        if(icon(1).eq.3.or.icon(1).eq.7) then
          nth0 = -ntheta
        else
          nth0 = ntheta
        endif
C
C       PLOT CROSS SECTION
C
	if(i .eq. 1) then
	  icon9 = icon(6)
	else
	  icon9 = icon(9)
	endif
        WRITE(ioutp,9901)'Sigma    ',I,LTRT(I)*JTRT(I),ALPHA
     1	,(IDAT(M),M=1,3)
        IF(I.eq.0 .and. eta(1).ne.0.0) then
          CALL DWPLOT(NTH0,ICON9 ,Csig(1)     ,DTHETA,THETA1)
        else
          CALL DWPLOT(NTH0,ICON9 ,Sigma(1,1),DTHETA,THETA1)
        endif
      NTEMP=-10
      M1=-1
      DO 2010 M=1,NTHETA
      IF(ABS(Sigma(M,3)).GT.0.01) then
        M1=0
        go to 2011
      ENDIF
 2010 CONTINUE
 2011 IF(IS(1).NE.0.AND.M1.EQ.0) THEN
c
c     Plot Asymmetries
c
        WRITE(ioutp,9901)'Asymmetry ',I,LTRT(I)*JTRT(I),ALPHA
     1	,(IDAT(M),M=1,3)
        CALL DWPLOT(NTH0,NTEMP  ,Sigma(1,3),DTHETA,THETA1)
      ENDIF
      ENDIF
 3000 CONTINUE
 3001 CONTINUE
      rewind 3
      RETURN
c
 9027 FORMAT('0  Theta    Sig(1)/Coul     Sigma(1)  ','    Pol   '
     1       ,'    Asy   ','    Ayy   ','    A22   ','    A21   '
     2       ,'    A20   ','   Theta  ','Cos(Theta)')
 9028 FORMAT('0  Theta      Sigma('     ,i2,   ')  ','    Pol   '
     1       ,'    Asy   ','    Ayy   ','    A22   ','    A21   '
     2       ,'    A20   ','   Theta  ','Cos(Theta)')
 9030 FORMAT(0PF8.2,2X,1P2E13.4,0P6F10.4,0PF10.2,0pf10.3)
 9031 FORMAT(0PF8.2,2X,1P1E13.4,0P6F10.4,0PF10.2,0pf10.3)
 9901 FORMAT('1',a12,' for channel',I2,10H   2*spin=,I3
     1,14X,15A4,I4,2(1H/,I2.2),I4,2(1H:,I2.2))
 9905 FORMAT(a9,1PE14.4)
      END

