C
C
C --------------------------------------------------------------
C
      SUBROUTINE Setr (Nent, Ntot, Ishift, Lpent, Lspin, Bound, Echan,
     *  Parext, Iflext, Beta, Sinsqr, Sin2ph, Dphi, Cscs, Sinphi,
     *  Cosphi, Zkfe, Zkte, Zeta, Alphar, Alphai, Not, Dpdr, Dsdr, Rmat,
     *  Ymat, Rootp, Elinvr, Elinvi, Psmall, Krext, Lrmat, Min, Max,
     *  Ipoten)
C
C *** PURPOSE -- GENERATE Linv = 1/(S-B+IP)
C ***                    Rootp = sqrt(P)
C ***                     Rmat = SUM Beta*Beta/((DEL E)-i(GAMGAM/2))
C ***                     Ymat = Linv - Rmat
C ***            Also return Lrmat = 1 if no R-matrix contribution
C
C       ROS 1999 update for Coulomb
C       NML most recent update September 2000
C
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      INCLUDE 'B1ZYX'
      INCLUDE 'B2ZYX'
      INCLUDE 'B3ZYX'
      INCLUDE 'B4ZYX'
      INCLUDE 'B44ZYX'
C
      DIMENSION Ishift(*), Lpent(*), Lspin(*), Bound(*), Echan(*),
     *   Parext(Krext,Ntotc,*), Iflext(Krext,Ntotc,*), Beta(Ntriag,*),
     *   Sinsqr(*), Sin2ph(*), Dphi(*), Cscs(2,*), Zkfe(*),
     *   Zkte(*), Zeta(*), Alphar(*), Alphai(*), Not(*), Dpdr(*),
     *   Dsdr(*), Rmat(2,*), Ymat(2,*), Rootp(*), Elinvr(*),
     *   Elinvi(*), Psmall(*)
      DIMENSION Sinphi(*), Cosphi(*)
c
C      DIMENSION Ishift(Ntotc), Lpent(Ntotc), Lspin(Ntotc), Bound(Ntotc),
C     *   Echan(Ntotc), Parext(Nrext,Ntotc,Ngroup),
C     *   Iflext(Nrext,Ntotc,Ngroup), Beta(Ntriag,Nres),
C     *   Zkte(Ntotc), Zeta(Ntotc), Alphar(Nres), Alphai(Nres),
C     *   Not(Nres), Dpdr(Ntotc), Dsdr(Ntotc), Rmat(2,Ntotc),
C     *   Ymat(2,Ntotc), Rootp(Ntotc), Elinvr(Ntotc), Elinvi(Ntotc),
C     *   Psmall(Ntotc)
C
      DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/, Tiny/1.d-8/
C
      IF (Nent.GT.10) STOP 'Nent > 10 in Setr in mxct3.f'
      Nnntot = Ntot + 1
      DO I=1,Ntot
         Nnntot = Nnntot - 1
         IF (Su.GE.Zero .AND. Su.LE.Echan(Nnntot)) THEN
            Ntot = Nnntot - 1
         ELSE
            GO TO 6
         END IF
      END DO
    6 CONTINUE
C
C *** INITIALIZE Rmat = R-MATRIX
C
      KL = 0
      DO K=1,Ntot
         IF (Numext.NE.0 .AND. Iflext(1,K,Nnnn).NE.-1) Aloge
     *      = dLOG( (Parext(2,K,Nnnn)-Su)/(Su-Parext(1,K,Nnnn)) )
         DO L=1,K
            KL = KL + 1
            Rmat(1,KL) = Zero
            Rmat(2,KL) = Zero
            IF (L.EQ.K .AND. Numext.NE.0 .AND. Iflext(1,K,Nnnn).NE.-1)
     *         THEN
               Rmat(1,KL) =   Parext(3,K,Nnnn)
     *                      + Parext(4,K,Nnnn)*SU
     *                      - Parext(5,K,Nnnn)*Aloge
               IF (Nrext.EQ.7) Rmat(1,KL) = Rmat(1,KL)
     *            + Parext(7,K,Nnnn)*Su**2 - Parext(6,K,Nnnn)*
     *                  (  Parext(2,K,Nnnn) - Parext(1,K,Nnnn)  )
     *            - Parext(6,K,Nnnn)*Aloge*Su
            END IF
         END DO
      END DO
C
      IF (Max.GE.Min .AND. Min.GT.0) THEN
         DO Ires=Min,Max
            KL = 0
            DO K=1,Ntot
               DO L=1,K
                  KL = KL + 1
                  IF (Su.GT.Echan(K) .AND. Su.GT.Echan(L) .AND.
     *              Beta(KL,Ires).NE.Zero) THEN
                    Rmat(1,KL) = Rmat(1,KL) + Alphar(Ires)*Beta(KL,Ires)
                    IF (Not(Ires).NE.1) THEN
                       Rmat(2,KL) = Rmat(2,KL) +
     *                    Alphai(Ires)*Beta(KL,Ires)
                    END IF
                  END IF
               END DO
            END DO
         END DO
      END IF
C
C *** Check if Rmat is Zero; if so, set Lrmat=1
      KL = 0
      DO K=1,Ntot
         DO L=1,K
            KL = KL + 1
            IF (Rmat(1,KL).NE.Zero) GO TO 90
            IF (Rmat(2,KL).NE.Zero) GO TO 90
         END DO
      END DO
      Lrmat = 1
   90 CONTINUE
C
C
C
C *** GENERATE Rootp,PH,(H+Rmat) MATRICES
C ***          Rootp = SQRT(P)
C ***             PH = 1/(S-B+IP)
C ***           Ymat = (1/(S-B+IP) - Rmat)
C ***           Dpdr = partial P wrt Rho
C ***           Dsdr = partial S wrt Rho
C
      KL = 0
      DO K=1,Ntot
         DO L=1,K
            KL = KL + 1
            Ymat(1,KL) = - Rmat(1,KL)
            Ymat(2,KL) = - Rmat(2,KL)
         END DO
      END DO
C
      CALL Zero_Array (Psmall, Ntot)
c
      KK = 0
      II = 0
      DO I=1,Ntot
         II = II + I
         Rootp (I) =  One
         Elinvr(I) =  Zero
         Elinvi(I) = -One
         Dpdr(I) = Zero
         Dsdr(I) = Zero
         Iffy = 0
         IF (SU.GT.Echan(I)) THEN
            IF (Lpent(I).NE.1) THEN
C
C ***          HERE PENETRABILITY IS Not CALCULATED
               Ymat(2,Ii) = Ymat(2,Ii) - One
C
            ELSE
C
               Lsp = Lspin(I)
               Q = dSQRT(Su-Echan(I))
               Rho = Zkte(I)*Q
C
               IF (Zeta(I).EQ.Zero) THEN
C ***             Here non-Coulomb phases and penetrabilities are 
C ***                to be calculated
                  IIII = I
                  CALL Sinsin (Lsp, Sinsqr(I), Sin2ph(I), Dphi(I),
     *               Zkfe(I), Sinphi(I), Cosphi(I))
                  CALL Pgh (Rho, Lsp, Bound(I), Hr, Hi, P, Dp, Ds,
     *               Ishift(I), Iffy)
C ***             HR and HI are real and imag parts of 1/(S-B+IP)
C ***                except when S-B+iP=Zero, in which case iffy=1
C
               ELSE IF (Zeta(I).NE.Zero) THEN
C ***             Here Coulomb penetrability is to be calculated
                  Eta  = Zeta(I)/Q
                  Jdopha = 1
                  CALL Pghcou (Rho, Lsp, Bound(I), HR, HI, P, DP, DS,
     *               Ishift(I), Iffy, Eta, Sinphi(i), Cosphi(i),
     *               Dphi(I), Jdopha)
                  Sinsqr(I) = Sinphi(I)**2
C                           = sin^2 (phase shift)
                  Sin2ph(I) = Two*Sinphi(I)*Cosphi(I)
C                           = sin(2 * phase shift)
               END IF
               IF (Ifdif.NE.0) THEN
                  Kk = Kk + I - 1
                  IF (I.GT.1) THEN
                     DO Kx=1,I-1
                        Kkx = Kk + Kx
                        Cscs(1,Kkx) = Cosphi(Kx)*Cosphi(I) - 
     *                                Sinphi(Kx)*Sinphi(I)
                        Cscs(2,Kkx) = Cosphi(Kx)*Sinphi(I) +
     *                                Sinphi(Kx)*Cosphi(I)
                     END DO
                  END IF
               END IF
C
cx               IF (Lrmat.NE.1) THEN
                  Rootp(I) = dSQRT(P)
                  Dpdr(I) = Dp
                  Dsdr(I) = Ds
                  IF (Iffy.EQ.0 .AND. .NOT. (Ishift(I).EQ.0 .AND.
     *               (One-P*Rmat(2,Ii).EQ.One .OR. P.LT.Tiny))) THEN
                     Elinvr(I) = Hr
                     Elinvi(I) = Hi
                     Ymat(1,Ii) = Hr + Ymat(1,Ii)
                     Ymat(2,Ii) = Hi + Ymat(2,Ii)
                  ELSE
c ***                here penetrability is very small but non-Zero
                     Psmall(I) = Rootp(I)
                     Ymat(1,Ii) = P*Ymat(1,Ii)
                     Ymat(2,Ii) = P*Ymat(2,Ii) - One
                     Rmat(1,Ii) = P*Rmat(1,Ii)
                     Rmat(2,Ii) = P*Rmat(2,Ii)
                     IF (Ntot.GT.1) THEN
                        IF (I.GT.1) THEN
                           DO J=1,I-1
                              Ji = (I*(I-1))/2 + J
                              Ymat(1,Ji) = Rootp(I)*Ymat(1,Ji)
                              Ymat(2,Ji) = Rootp(I)*Ymat(2,Ji)
                              Rmat(1,Ji) = Rootp(I)*Rmat(1,Ji)
                              Rmat(2,Ji) = Rootp(I)*Rmat(2,Ji)
                           END DO
                        END IF
                        IF (I.LT.Ntot) THEN
                           DO J=I+1,Ntot
                              Ji = (J*(J-1))/2 + I
                              Ymat(1,Ji) = Rootp(i)*Ymat(1,Ji)
                              Ymat(2,Ji) = Rootp(i)*Ymat(2,Ji)
                              Rmat(1,Ji) = Rootp(i)*Rmat(1,Ji)
                              Rmat(2,Ji) = Rootp(i)*Rmat(2,Ji)
                           END DO
                        END IF
                     END IF
                     Rootp(I) = One
                     Elinvr(I) = Zero
                     Elinvi(I) = -One
                     Dpdr(I) = Zero
                     Dsdr(I) = Zero
                  END IF
cx               END IF
            END IF
         ELSE
            Rootp(I)  = Zero
            Elinvr(I) = One
            Elinvi(I) = Zero
            Dpdr(I) = Zero
            Dsdr(I) = Zero
         END IF
      END DO
C
      IF (Lrmat.NE.1) THEN
C ***    Check if one channel is irrelevant; if so, set to unity
         KL = 0
         DO K=1,Ntot
            KL = KL + K
            IF (Ymat(1,KL).EQ.Zero .AND. Ymat(2,KL).EQ.Zero) THEN
               Ymat(1,KL) = One
            END IF
         END DO
      END IF
      RETURN
      END
C
C
C --------------------------------------------------------------
C
      SUBROUTINE Setxqx (Ntot, Yinv, Rmat, Xqr, Xqi, Rootp, Elinvr,
     *   Elinvi, Xxxxr, Xxxxi)
C
C *** purpose -- form XQ & XXXX matrices, where
C ***            XQ   = Yinv * Rmat       and
C ***            XXXX = P/L + sqrt(P)/L         (1/L-R)**-1          sqrt(P)/L
C ***                 =       sqrt(P)/(S-B+IP) * Yinv       * Rmat * sqrt(P)
C ***                 =       sqrt(P)/L        * XQ                * sqrt(P)
C
C ***            Note that the matrix W defined in SAMMY manual is given
C ***                 by W(c,c') = delta(c,c') + 2i XXXX(c,c')
C ***                 as in Eq. (III.D.4) in SAMMY manual R3
C
C ***         ie W    = I + 2i XXXX
C
C     NML, January 2000
C
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      INCLUDE 'B1ZYX'
      INCLUDE 'B2ZYX'
      INCLUDE 'B44ZYX'
C
      DIMENSION Yinv(2,*), Rmat(2,*), Xqr(Ntot,*), Xqi(Ntot,*),
     *   Rootp(*), Elinvr(*), Elinvi(*), Xxxxr(*), Xxxxi(*)
C      DIMENSION Yinv(2,Nnnn), Rmat(2,Nnnn), Xqr(Ntot,Ntot),
C     *   Xqi(Ntot,Ntot), Rootp(Ntotc), Elinvr(Ntot), Elinvi(Ntot),
C     *   Xxxxr(Nnnn), Xxxxi(Nnnn)
      EXTERNAL ijkl
C
      CALL Zero_Array (Xqr, Ntot*Ntot)
      CALL Zero_Array (Xqi, Ntot*Ntot)
C
C *** Xqr(k,i) = (L**-1-R)**-1 * R ... note asymmetry
      DO I=1,Ntot
         DO J=1,Ntot
            Ij = Ijkl(J,I)
            DO K=1,Ntot
               Jk = Ijkl(K,J)
               Xqr(K,I) = Xqr(K,I) + Yinv(1,Ij)*Rmat(1,Jk) -
     *                               Yinv(2,Ij)*Rmat(2,Jk)
               Xqi(K,I) = Xqi(K,I) + Yinv(1,Ij)*Rmat(2,Jk) +
     *                               Yinv(2,Ij)*Rmat(1,Jk)
            END DO
         END DO
      END DO
C
C *** Xxxx = sqrt(P)/L  * xq * sqrt(P) ... symmetric
      IJ = 0
      DO I=1,Ntot
         Plr = Rootp(I)*Elinvr(I)
         Pli = Rootp(I)*Elinvi(I)
         DO J=1,I
            Ij = Ij + 1
            Xxxxr(Ij) = Rootp(J)* (Xqr(J,I)*Plr-Xqi(J,I)*Pli)
            Xxxxi(Ij) = Rootp(J)* (Xqi(J,I)*Plr+Xqr(J,I)*Pli)
         END DO
      END DO
      RETURN
      END
C
C
C --------------------------------------------------------------
C
      SUBROUTINE Sectio (Nent, Next, Lspin, If_Excl, Ifcros, Zke, Zeta,
     *   Xxxxr, Xxxxi, Sinsqr, Sin2ph, Termf, Crss, Crssx, Cscs, Dgoj,
     *   Ntotnn)
C
C     NML, April 2000
C *** purpose -- generate pieces of cross sections (except for "4 pi / E")
C
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      INCLUDE 'B1ZYX'
      INCLUDE 'B2ZYX'
      INCLUDE 'B44ZYX'
C
      DIMENSION If_Excl(*), Ifcros(*), Zke(*), Xxxxr(*), Xxxxi(*),
     *   Sinsqr(*), Sin2ph(*), Termf(*), Crss(*),
     *   Crssx(2,Ntotc,Ntotc,*), Cscs(2,*), Lspin(*), Zeta(*)
      DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
C
C
C
      IF (Ifcros(1).EQ.1) THEN
C ***    elastic  Crss(1) = g*0.25* sum(entrance chs c,c')
C ***                                     times |(1-U(c,c'))| **2 / Zz
C ***                     = g* [ sin(phi)**2 * (1-2XXXXi)
C ***                            - sin(2phi)*XXXXr
C ***                            + (XXXXr**2 + XXXXi**2) ] / Zz
         Crss(1) = Zero
         Ii = 0
         Ij = 0
         DO I=1,Nent
            Zz = Zke(I)**2
            Ii = Ii + I
            Termn = Sinsqr(I)*( One - Two * Xxxxi(Ii) )
     *                    - Sin2ph(I)*Xxxxr(Ii)
            Termn = Termn / Zz
            DO J=1,I
               Ij = Ij + 1
               Ar = ( Xxxxr(Ij)**2 + Xxxxi(Ij)**2 )/Zz
               IF (I.NE.J) Ar = Ar + Ar
               Termn = Termn + Ar
            END DO
            Crss(1) = Termn + Crss(1)
         END DO
         Crss(1) = Crss(1)*Dgoj
      END IF
C
      IF (Ifcros(2).EQ.1) THEN
C ***    absorption = g*0.25 * sum(inc c) 
C ***                  [ 1 -  sum(inc c') |U(c,c')| **2 ] / Zz
C ***               = - g* (XXXXr**2 + XXXXi**2) / Zz
         Crss(2) = Zero
         Ii = 0
         Ij = 0
         DO I=1,Nent
            Ii = Ii + I
            Zz = Zke(I)**2
            Terma =  Xxxxi(Ii) / Zz
            DO J=1,I
               Ij = Ij + 1
               Ar = (- Xxxxr(Ij)**2 - Xxxxi(Ij)**2) / Zz
               IF (I.NE.J) Ar = Ar + Ar
               Terma = Terma + Ar
            END DO
            Crss(2) = Terma + Crss(2)
         END DO
         Crss(2) = Crss(2)*Dgoj
      END IF
C
      IF (Next.GT.0 .and. Ntotnn.GT.Nent) THEN
C ***    reaction ch c'= g*0.25 * sum(inc c) |U(c,c')|**2 / Zz
C ***                  = g* (XXXXr**2 + XXXXi**2) / Zz
         DO Jj=1,Next
            Termf(Jj) = Zero
         END DO
         DO Jj=1,Next
            IF (Jj+Nent.LE.Ntotnn .AND.
     *         Ifcros(Jj+2).NE.0 .AND. If_Excl(jj+Nent).EQ.0) THEN
               J = Jj + Nent
               DO I=1,Nent
                  Zz = Zke(I)**2
                  Ij = (J*(J-1))/2 + I
Cq                Ij = Ijkl(i,j) but i < j always
                  Termf(Jj) = Termf(Jj) + (Xxxxr(Ij)**2+Xxxxi(Ij)**2)/Zz
               END DO
               Crss(2+Jj) = Termf(Jj)*Dgoj
            END IF
         END DO
      END IF
C
C
      IF (Ifdif.NE.0) THEN
C
C ***    Angular Distribution   Crssx(.,i,ix,igroup) = (1-U)/2
         DO Ichan=1,Ntotnn
            IF (If_Stay (Ichan, Ifdif, Nent, If_Excl(Ichan)) .EQ.0) THEN
C
               IF (Zeta(Ichan).NE.Zero) THEN
                  CALL Get_Coul_Phase (Cr, Ci, Lspin(Ichan),
     *               Zeta(Ichan), Su)
               ELSE
                  Cr = One
                  Ci = Zero
               END IF
               II = (Ichan*(Ichan-1))/2
               DO Ichanx=1,Ichan
                  II = II + 1
                  IF (Ichanx.LE.Nent) THEN
C ***                real and imaginary parts of (1-U)/2
                     IF (Ichanx.EQ.Ichan) THEN
                        Ar = Sinsqr(Ichan)*(One-Two*Xxxxi(II))
     *                          - Sin2ph(Ichan)*Xxxxr(II) + Xxxxi(II)
                        Ai = Sin2ph(Ichan)*(0.5d0-Xxxxi(II))
     *                          - (One-Two*Sinsqr(Ichan)) * Xxxxr(II)
                     ELSE
                        Ar = Cscs(1,II)*Xxxxi(II) - Cscs(2,II)*Xxxxr(II)
                        Ai =-Cscs(1,II)*Xxxxr(II) - Cscs(2,II)*Xxxxi(II)
                     END IF
                     If (Zeta(Ichan ).NE.Zero .OR.
     *                   Zeta(Ichanx).NE.Zero) THEN
                        Br = Ar*Cr - Ai*Ci
                        Bi = Ar*Ci + Ai*Cr
                        IF ((Lspin(Ichanx).NE.Lspin(Ichan) .OR.
     *                     Zeta(Ichanx).NE.Zeta(Ichan)) .AND.
     *                     Ichan.NE.Ichanx ) THEN
                           IF (Zeta(Ichanx).NE.Zero) THEN
                              CALL Get_Coul_Phase (Dr, Di,
     *                           Lspin(Ichanx), Zeta(Ichanx), Su)
                           ELSE
                              Dr = One
                              Di = Zero
                           END IF
                        ELSE
                           Dr = Cr
                           Di = Ci
                        END IF
                        Ar = Br*Dr - Bi*Di
                        Ai = Br*Di + Bi*Dr
                     END IF
                     Crssx(1,Ichanx,Ichan,Nnnn) = Ar
                     Crssx(2,Ichanx,Ichan,Nnnn) = Ai
                  END IF
               END DO
            END IF
         END DO
      END IF
C
C
      RETURN
      END
C
C
C --------------------------------------------------------------
C
      INTEGER FUNCTION Ijkl (M,N)
      IF (M.LE.N) THEN
         Ijkl = (N*(N-1))/2 + M
      ELSE
         Ijkl = (M*(M-1))/2 + N
      END IF
      RETURN
      END
C
C
C --------------------------------------------------------------
C
      SUBROUTINE Get_Coul_Phase (Cr, Ci, Lspin, Zeta, Su)
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      DATA Zero /0.0d0/, One /1.0d0/
      Eta = Zeta/Dsqrt(SU)
      L = Lspin
      IF (L.EQ.0) THEN
         Cr = One
         Ci = Zero
      ELSE
         aa = eta**2 + One
         aa = dsqrt(aa)
         Cc = One/aa
         Ss = eta/aa
         IF (L.GT.1) THEN
            DO LL=2,L
               aa = eta**2 + LL**2
               aa = dsqrt(aa)
               ccx = Ll/aa
               ssx = eta/aa
               ccy = cc*ccx - ss*ssx
               ssy = cc*ssx + ss*ccx
               cc = ccy
               ss = ssy
            END DO
         END IF
         Cr = Cc
         Ci = Ss
C ***    (Cr,Ci) = (Re,Im) exp{ 
C ***                    i sum from (LL=1 to L) tan^{-1} (eta/LL) }
      END IF
      RETURN
      END
C
C
C --------------------------------------------------------------
C
      INTEGER FUNCTION If_Stay (Ichan, Ifdif, Nent, If_Excl)
C *** If_Stay = 0 if want this channel, If_Stay = 1 if do not
      If_Stay = 0
      IF (Ifdif.EQ.1 .AND. Ichan.GT.Nent) THEN
         If_Stay = 1
      ELSE IF (Ifdif.EQ.2) THEN
         IF (Ichan.LE.Nent) THEN
            If_Stay = 1
         ELSE IF (If_Excl.EQ.1) THEN
            If_Stay = 1
         END IF
      END IF
      RETURN
      END