-
Since this is file move only, this commit will not compile
Since this is file move only, this commit will not compile
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
mxct09.f90 9.33 KiB
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. (II.B1.3) in SAMMY manual R7
C
C *** ie W = I + 2i XXXX
C
use fixedi_m
use ifwrit_m
use varyr_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
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, igr, Echan, If_Excl, Ifcros, Zke,
* Zeta, Xxxxr, Xxxxi, Sinsqr, Sin2ph, Termf, Crss, Crssx, Cscs,
* Dgoj, Ntotnn)
C
C *** Purpose -- Generate pieces of cross sections (except for "4 pi / E")C
use fixedi_m
use ifwrit_m
use varyr_common_m
use EndfData_common_m
use SammyResonanceInfo_M
use RMatResonanceParam_M
IMPLICIT DOUBLE PRECISION (a-h,o-z)
C
type(SammySpinGroupInfo)::spinInfo
type(SammyChannelInfo)::channelInfo
type(RMatChannelParams)::channel
DIMENSION Echan(*), If_Excl(*), Ifcros(*), Zke(*),
* Zeta(*), Xxxxr(*), Xxxxi(*),Sinsqr(*), Sin2ph(*), Termf(*),
* Crss(Ncrsss), Crssx(2,Ntotc,Ntotc,*), Cscs(2,*)
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)*(One-Xxxxi(Ii)) - Xxxxr(Ii)**2
DO J=1,I
Ij = Ij + 1
IF (J.NE.I) THEN
Ar = - Xxxxr(Ij)**2 - Xxxxi(Ij)**2
Ar = Ar + Ar
Terma = Terma + Ar
END IF
END DO
Terma = Terma / Zz
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) THEN
IF ( (Kaptur.EQ.0 .AND. If_Excl(Jj+Nent).EQ.0) .OR.
* (Kaptur.EQ.1 .AND. If_Excl(Jj+Nent).EQ.1) ) 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
END IF
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
call resParData%getSpinGroupInfo(spinInfo, igr)
DO Ichan=1,Ntotnn
call spinInfo%getChannelInfo(channelInfo, Ichan)
call resParData%getChannel(channel, channelInfo)
Lspin = channel%getL()
Ifs = If_Stay (Ichan, Ifdif, Nent, If_Excl(Ichan), Kaptur)
IF (Ifs.EQ.0) THEN
C
IF (Zeta(Ichan).NE.Zero .AND. Su.GT.Echan(Ichan)) THEN
CALL Get_Coul_Phase (Cr, Ci, Lspin,
* Echan(Ichan), Zeta(Ichan), Su)
ELSE
Cr = One
Ci = Zero
END IF
DO Ichanx=1,Nent
call spinInfo%getChannelInfo(channelInfo, Ichanx)
call resParData%getChannel(channel, channelInfo)
Lspinx = channel%getL()
IF (Ichanx.LE.Ichan) THEN
II = (Ichan*(Ichan-1))/2 + Ichanx
ELSE
II = (Ichanx*(Ichanx-1))/2 + Ichan
END IF
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 ((Lspinx.NE.Lspin .OR.
* Zeta(Ichanx).NE.Zeta(Ichan)) .AND.
* Ichan.NE.Ichanx ) THEN
IF (Zeta(Ichanx).NE.Zero .AND. * Su.GT.Echan(Ichanx)) THEN
CALL Get_Coul_Phase (Dr, Di,
* Lspinx, Echan(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
C entrance,exit ,group
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
INTEGER FUNCTION If_Stay (Ichan, Ifdif, Nent, If_Excl, Kaptur)
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
C *** Ifdif=1 means elastic
If_Stay = 1
ELSE IF (Ifdif.EQ.2) THEN
C *** Ifdif=2 means reaction of some kind
IF (Ichan.LE.Nent) THEN
If_Stay = 1
C *** Do not want elastic
ELSE IF (Kaptur.EQ.0 .AND. If_Excl.EQ.1) THEN
If_Stay = 1
C *** Do not want excluded channel in final state
ELSE IF (Kaptur.EQ.1 .AND. If_Excl.EQ.0) THEN
If_Stay = 1
C *** Will subtract only excluded channels from absorption
ELSE IF (If_Excl.EQ.-1) THEN
If_Stay = -1
C *** Do not want excluded channel anywhere in the calculation
END IF
END IF
RETURN
END