From 19f6d2bcbc07fefda7ec910fb16586d9d818d706 Mon Sep 17 00:00:00 2001
From: Wiarda <wiardada@ornl.gov>
Date: Fri, 10 Sep 2021 12:00:37 -0400
Subject: [PATCH] Convert MLBW to use the new frame work
---
sammy/src/cro/mcro2a.f | 5 +-
sammy/src/mlb/mmlb0.f | 166 ----------
sammy/src/mlb/mmlb0.f90 | 59 ++++
sammy/src/mlb/mmlb1.f | 165 ----------
sammy/src/mlb/mmlb2.f | 32 --
sammy/src/mlb/mmlb2.f90 | 39 +++
sammy/src/mlb/mmlb3.f | 227 -------------
sammy/src/mlb/mmlb3.f90 | 242 ++++++++++++++
sammy/src/mlb/mmlb4.f | 352 ---------------------
sammy/src/mlb/mmlb4.f90 | 287 +++++++++++++++++
sammy/src/sammy/CMakeLists.txt | 9 +-
sammy/src/the/CrossSectionCalcDriver_M.f90 | 2 +
sammy/src/the/DerivativeHandler.cpp | 5 +-
sammy/src/the/mthe0.f90 | 1 +
sammy/src/xxx/mxxx6.f90 | 15 +-
15 files changed, 648 insertions(+), 958 deletions(-)
delete mode 100644 sammy/src/mlb/mmlb0.f
create mode 100644 sammy/src/mlb/mmlb0.f90
delete mode 100644 sammy/src/mlb/mmlb1.f
delete mode 100755 sammy/src/mlb/mmlb2.f
create mode 100755 sammy/src/mlb/mmlb2.f90
delete mode 100644 sammy/src/mlb/mmlb3.f
create mode 100644 sammy/src/mlb/mmlb3.f90
delete mode 100644 sammy/src/mlb/mmlb4.f
create mode 100644 sammy/src/mlb/mmlb4.f90
diff --git a/sammy/src/cro/mcro2a.f b/sammy/src/cro/mcro2a.f
index 18e89f759..b660be6c8 100644
--- a/sammy/src/cro/mcro2a.f
+++ b/sammy/src/cro/mcro2a.f
@@ -253,8 +253,9 @@ C
IF (Kcros.LE.2) THEN
C *** CALCULATE SIN AND COS OF POTENTIAL SCATTERING PHASE SHIFT,
C *** AND R-EXTERNAL PHASE SHIFT
- CALL Cossin (Zke(1,N),
- * A_Ics, A_Isi, A_Idphi, Nnnn, Ipoten)
+ CALL Cossin (resparData, Zke(1,N),
+ * A_Ics, A_Isi, A_Idphi, Nnnn, Ipoten,
+ * Squ, Su)
END IF
C
C
diff --git a/sammy/src/mlb/mmlb0.f b/sammy/src/mlb/mmlb0.f
deleted file mode 100644
index 624e48073..000000000
--- a/sammy/src/mlb/mmlb0.f
+++ /dev/null
@@ -1,166 +0,0 @@
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Sammlb_0
-C
-C *** CAUTION *** Sammlb has not been updated as well as should perhaps be,
-C *** which may not matter very much because its use is not
-C *** really recommended anyway. Users are advised to proceed
-C *** with caution when using any "new" features with MLB.
-C
- use oops_common_m
- use fixedi_m
- use ifwrit_m
- use exploc_common_m
- use array_sizes_common_m
- use oopsch_common_m
- use cbro_common_m
- use lbro_common_m
- use templc_common_m
- use AllocateFunctions_m
- use rsl7_m
- use xct_m
- use mxct27_m
- use EndfData_common_m, only : covData, resParData, radFitFlags
- use AuxGridHelper_M, only : setAuxGridOffset, setAuxGridRowMax
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- real(kind=8),allocatable,dimension(:)::A_Idum
-C
-C
-C
- CALL Initil
- IF (Krmatx.EQ. 1) WRITE (6,99999)
-99999 FORMAT (' *** SAMMY-MLBW 15 Nov 07 *** multilevel Breit Wigner')
- IF (Krmatx.EQ.-1) WRITE (6,99998)
-99998 FORMAT (' *** SAMMY-MLBW 15 Nov 07 ** single level Breit Wigner')
- WRITE (21,99997)
-99997 FORMAT (' *** Caution *** This version of SAMMY-MLB has been teste
- *d to be sure',
- * /, ' cross sections are correct, but no guara
- *ntees are made',
- * /, ' about the derivatives wrt radii or abund
- *ances, e.g.',
- * /, ' See NML if you suspect trouble here.')
-C
- Segmen(1) = 'M'
- Segmen(2) = 'L'
- Segmen(3) = 'B'
- Nowwww = 0
-C
- IF (covData%getPupedParam().GT.0) THEN
- WRITE (6,10100)
-10100 FORMAT ('Breit-Wigner coding does not include options for', /,
- * 'PUPs (Propagated-Uncertainty Parameters, Flag=3).', /,
- * 'Contact N.M.Larson at LarsonNM@ornl.gov if you', /,
- * 'desperately need this option.')
- END IF
-C
-C *** organize for broadening
-C
- Icr = Kcros
- Npfil3 = 1
- Nrfil3 = 1
-C
-C *** guesstimate size of array needed for SAMMY-MLBW
- call setAuxGridOffset(1) ! reset starting point for auxillary grid
- call setAuxGridRowMax(0)
- CALL Estmlb (N1, Ntwo, N2, N3, N5, N6, N7, Nrfil3, Npfil3)
-C
-C ### one ###
- Ks_Res = Ksolve
- CALL Set_Kws_Xct
- call calcData%setUpList(resParData, radFitFlags, Iq_Iso)
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - <
-C
-C ### three ###
- N = Ndatb
- N = Ntotc*Ngroup
- call allocate_real_data(A_Iaaone, N)
- call allocate_real_data(A_Iaatwo, N)
- N = Ntriag*Ngroup
- call allocate_real_data(A_Iaathr, N)
- N = Ntotc
- call allocate_real_data(A_Ics, N) ! this was set equal to Ip
- call allocate_real_data(A_Isi, N) ! this was set equal to Igami
- call allocate_real_data(A_Idphi, N)
- N = N6
- call allocate_real_data(A_Ipiece, N)
- call allocate_real_data(A_Idum, N7)
-C
-C
-C
-C *** Sub Work generates Theory and derivatives
- CALL Work_Mlb ( I_Iflmsc , A_Iprnbk , I_Iflnbk , A_Iprbgf ,
- * I_Iflbgf , I_Indbgf , A_Ibgfmi , A_Ibgfma , A_Itexbg ,
- * A_Iteabg , A_Ith ,
- * A_Isigxx , A_Idasig , A_Idbsig ,
- * A_Ipiece , A_Idum , A_Iadder , A_Iaddcr , I_Inbt ,
- * I_Iint)
- deallocate(A_Ics)
- deallocate(A_Isi)
- deallocate(A_Idphi)
- deallocate(A_Iaaone)
- deallocate(A_Iaatwo)
- deallocate(A_Iaathr)
- deallocate(A_Idum)
-C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ->
-C
-C
-C *** if making fake data write it out, then stop
- IF (Kfake.EQ.1) CALL Newdat(A_Ith )
- IF (Kfake.EQ.1) STOP '[STOP in Sammlb_0 in mlb/mmlb0.f]'
-C
-C
- Nnnsig = 1
- Kiniso = Kkkiso
- Niniso = Nnniso
- RETURN
-C
- END
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Estmlb (N1, Ntwo, N2, N3, N5, N6, N7, Nr, Np)
-C *** guesstimate size of array needed for SAMMY - Multi Level Breit Wigner
-C
- use fixedi_m
- use ifwrit_m
- use broad_common_m
- use lbro_common_m
- use rsl7_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
- N1 = Ntotc*Nfpres
- IF (N1.EQ.0) N1 = 1
-C ### one ###
- CALL Figure_Kws_Cro (Kone)
- K1 = Kone + Nres
-C
-C ### three ###
- Ntriag = (Ntotc*(Ntotc+1))/2
- N2 = 1
- N2 = Ntriag*Nfpres
- N3 = 1
- IF (Ksolve.NE.2) N3 = Ngbout
- N5 = 1
- IF (Numiso.NE.0) N5 = Numiso
- IF (N5.EQ.0) N5 = 1
- N6 = 1
- IF (Kpiece.EQ.1) N6 = Ndatb*Ngroup
- IF (N6.EQ.0) N6 = 1
- N7 = 1
- IF (Kpiece.EQ.1) N7 = Ndatb
- IF (N7.EQ.0) N7 = 1
- K = Ngroup + Ndatb + (Ntotc*2+Ntriag)*Ngroup + 2*N1 + N2
- * + N3 + 3*Ntotc + 3*N5 + N6 + N7
- IF (Jcros.EQ.6) K = K + N3
- K = K + K1
- k = k + 2*nr + 2*np
- K = Idimen (K, 1, 'K, 1')
- I = Idimen (K, -1, 'K, -1')
- I = Idimen (0, 0, '0, 0')
-C
- RETURN
- END
diff --git a/sammy/src/mlb/mmlb0.f90 b/sammy/src/mlb/mmlb0.f90
new file mode 100644
index 000000000..1775eba32
--- /dev/null
+++ b/sammy/src/mlb/mmlb0.f90
@@ -0,0 +1,59 @@
+module MlbwCalc_m
+ contains
+!
+!
+! --------------------------------------------------------------
+!
+ SUBROUTINE Sammlb_0
+!
+! *** CAUTION *** Sammlb has not been updated as well as should perhaps be,
+! *** which may not matter very much because its use is not
+! *** really recommended anyway. Users are advised to proceed
+! *** with caution when using any "new" features with MLB.
+!
+ use fixedi_m, only : Kiniso, Kkkiso, Niniso, Nnnsig, Nnniso
+ use ifwrit_m, only : Krmatx, Ks_Res, Ksolve
+ use oopsch_common_m, only : Nowwww, Segmen
+ use cbro_common_m, only : Segnam
+ use rsl7_m
+ use AuxGridHelper_M, only : setAuxGridOffset, setAuxGridRowMax
+ IMPLICIT None
+ integer::I,Kone, Idimen
+
+ external Idimen
+
+
+ IF (Krmatx.EQ. 1) WRITE (6,99999)
+99999 FORMAT (' *** SAMMY-MLBW 15 Nov 07 *** multilevel Breit Wigner')
+ IF (Krmatx.EQ.-1) WRITE (6,99998)
+99998 FORMAT (' *** SAMMY-MLBW 15 Nov 07 ** single level Breit Wigner')
+ WRITE (21,99997)
+99997 FORMAT (' *** Caution *** This version of SAMMY-MLB has been tested to be sure', &
+ /, ' cross sections are correct, but no guarantees are made', &
+ /, ' about the derivatives wrt radii or abundances, e.g.', &
+ /, ' See NML if you suspect trouble here.')
+!
+ Segmen(1) = 'M'
+ Segmen(2) = 'L'
+ Segmen(3) = 'B'
+ Nowwww = 0
+!
+! *** guesstimate size of array needed for SAMMY-MLBW
+ call setAuxGridOffset(1) ! reset starting point for auxillary grid
+ call setAuxGridRowMax(0)
+
+ CALL Figure_Kws_Cro (Kone)
+!
+! ### one ###
+ Ks_Res = Ksolve
+ CALL Set_Kws_Xct
+
+!
+!
+ Nnnsig = 1
+ Kiniso = Kkkiso
+ Niniso = Nnniso
+ RETURN
+!
+ END
+end module MlbwCalc_m
diff --git a/sammy/src/mlb/mmlb1.f b/sammy/src/mlb/mmlb1.f
deleted file mode 100644
index 50fe99457..000000000
--- a/sammy/src/mlb/mmlb1.f
+++ /dev/null
@@ -1,165 +0,0 @@
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Work_Mlb ( Iflmsc, Parnbk, Iflnbk, Parbgf,
- * Iflbgf, Kndbgf, Bgfmin, Bgfmax, Texbgf, Teabgf,
- * Theory, Sigxxx, Dasigx, Dbsigx,
- * Pieces, Dum, Adderg, Addcro, Nbt, Int)
-C
-C *** Purpose -- Generate theoretical data and partial derivatives
-C
- use oops_common_m
- use fixedi_m
- use ifwrit_m
- use samxxx_common_m
- use oopsch_common_m
- use fixedr_m
- use varyr_common_m
- use cbro_common_m
- use lbro_common_m
- use mxct27_m
- use SammyGridAccess_M
- use EndfData_common_m
- use array_sizes_common_m, only : calcData
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
- LOGICAL Ywhich
-C
- DIMENSION Iflmsc(*), Parnbk(*), Iflnbk(*), Parbgf(*), Iflbgf(*),
- * Kndbgf(*), Bgfmin(*), Bgfmax(*), Texbgf(*), Teabgf(*),
- * Theory(Nnnsig,*), Pieces(Ngroup,*),
- * Dum(*), Adderg(*), Addcro(*), Nbt(*), Int(*)
- DIMENSION Sigxxx(Nnnsig,*),
- * Dasigx(Nnnsig,*), Dbsigx(Nnnsig,Ndbxxx,*)
- type(SammyGridAccess)::grid
-C
-C DIMENSION Iflmsc(nummsc), Parnbk(Numnbk), Iflnbk(Numnbk),
-C * Theory(Nnnsig*Ndat), Pieces(Ngroup,Ndatb), Dum(Ndatb),
-C * Adderg(NP), Addcro(NP), Nbt(Nrfil3), Int(Nrfil3)
-C
- DATA Zero /0.0d0/
-C
- call grid%initialize()
- call grid%setParameters(numcro, ktzero)
- call grid%setToAuxGrid(expData)
- numEl = grid%getNumEnergies(expData)
- Ywhich = Ydoppr.OR.Yresol
-C
-C
- IF (Kpiece.EQ.1) CALL Zero_Array (Pieces, Ngroup*numEl)
- Kkkkkk = 0
- Kkkmin = 0
- call calcData%nullify()
- call calcData%setNnsig(Nnnsig)
-C
-C *** Calculate cross sections and derivatives for all
-C *** positive energies
- DO J=1,numEl
- Iskip = 0
- Jdat = J
- do Jsig = 1,Nnnsig
- Jcount = (Jdat-1)*Nnnsig + Jsig
- call calcData%setToZero(Jcount, Ndasig + Ndbsig)
- Theory(Jsig, Jdat) = 0.0d0
- end do
- CALL Zero_Array (Sigxxx, Nnnsig*Nnniso)
- IF (Ndasig.GT.0) CALL Zero_Array (DaSigx, Nnnsig*Ndasig)
- IF (Ndbsig.GT.0) CALL Zero_Array (DbSigx,
- * Nnnsig*Ndbsig*Nnniso)
- Su = grid%getEnergy(J, expData)
- IF (Kkkdop.NE.1 .OR. Su.GT.Zero) THEN
- Kkkkkk = Kkkkkk + 1
- IF (Su.NE.Zero) THEN
- IF (Su.GT.Emax .AND. Kartgd.EQ.1) THEN
- Sigxxx(1,1) = Zero
- GO TO 20
- END IF
-C
- IF (Su.LT.Zero) Su = - Su
- Squ = dSQRT(Su)
- Ipoten = 0
- IF (Kpoten.NE.0 .AND. (J.EQ.1 .OR.
- * J.EQ.numEl)) THEN
-10000 FORMAT(/, ' Potential Scattering for energy',
- * 1PE13.6, /,
- * ' group l for this group for this l')
- WRITE (6,10000) Su
- WRITE (21,10000) Su
- Ipoten = 1
- END IF
-C
-C ********* Generate cross sections and derivatives
- ee = grid%getEnergy(J, expData)
- IF (ee.NE.Zero) THEN
- CALL Zparsh_Mlb (Ipoten)
- IF (ee.LT.Zero) THEN
- CALL Fix_Negative_E (Sigxxx, Dasigx, Dbsigx)
- END IF
- END IF
- 20 CONTINUE
- END IF
-C
-C ********* If this is fake data, store it and move on
- IF (Kfake.EQ.1) THEN
- Theory(1,J) = Sigxxx(1,1)
- Iskip = 2
- END IF
-C
-C
-C
-C ********* If want Maxwellian averages or energy-averages with no
-C ********* correction terms, do same as if there is broadening
-C ********* except don't convert to transmission unless needed
- IF (Maxwel.EQ.1 .OR. Knocor.EQ.1) THEN
- Iskip = 1
- END IF
-C
-C ********* If there is Doppler broadening, skip the rest of this loop
- IF (Ydoppr) Iskip = 1
-C
-C ********* If no Doppler, and this is transmission, make conversion
-C
-C ********* If there is resolution, we're done here
- IF (Yresol) Iskip = 1
-C
-C
- IF (Iskip.EQ.0) THEN
-C
-C
-C ********* If ETA is being calculated and there is no broadening,
-C ********* find derivative with respect to nu
-C
-C ********* If there is normalization or background ...
-C
-C ********* Write results onto Theory if there is no broadening etc
- IF (Jjjdop.NE.1) THEN
- ! mlb doesn't have isotopes
- call calcData%addCalculatedData(Kkkkkk, Nnnsig,
- * ndasig, ndbsig, -1, Sigxxx,
- * Dasigx, Dbsigx)
- do Jsig = 1,Nnnsig
- Jcount = (Kkkkkk-1)*Nnnsig + Jsig
- Theory(Jsig,Jdat) =
- * calcData%getData(Jcount, 0, 1)
- end do
- END IF
-C
- ELSE IF (Iskip.EQ.1) THEN
-C ********* Here there is broadening to be done yet
-C
-C ********* Copy results
- IF (Kkkmin.EQ.0) Kkkmin = Jdat
- ! cro doesn't have isotopes
- call calcData%addCalculatedData(Jdat, Nnnsig, ndasig,
- * ndbsig, -1, Sigxxx, Dasigx, Dbsigx)
-C
- END IF
- END IF
- END DO
-C *** End of do-loop on energies
-C
-C
- IF (Kpiece.EQ.1) CALL Odfpcs (Pieces, Dum)
-C
- call grid%destroy()
- RETURN
- END
diff --git a/sammy/src/mlb/mmlb2.f b/sammy/src/mlb/mmlb2.f
deleted file mode 100755
index b0294bc10..000000000
--- a/sammy/src/mlb/mmlb2.f
+++ /dev/null
@@ -1,32 +0,0 @@
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Zparsh_Mlb (Ipoten)
-C
-C *** Purpose -- Form the cross section Sigma and the
-C *** ( Partial derivatives of the cross section with
-C *** respect to the varied parameters ) = Sjkl
-C
- use over_common_m
- use oops_common_m
- use fixedi_m
- use ifwrit_m
- use exploc_common_m
- use array_sizes_common_m
- use templc_common_m
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
-C
- CALL Abpart_Mlb (A_Iechan ,
- * A_Izke , A_Iaaone , A_Iaatwo, A_Iaathr ,
- * A_Ics , A_Isi , A_Idphi )
-C
-C
- CALL Parsh_Mlb (
- * A_Izke ,
- * A_Iaaone , A_Iaatwo , A_Iaathr , A_Ipiece , Ipoten,
- * A_Isigxx , A_Idasig , A_Idbsig)
-C
- RETURN
- END
diff --git a/sammy/src/mlb/mmlb2.f90 b/sammy/src/mlb/mmlb2.f90
new file mode 100755
index 000000000..de615e707
--- /dev/null
+++ b/sammy/src/mlb/mmlb2.f90
@@ -0,0 +1,39 @@
+module mmlb2_m
+ use MlbwCrossCalc_M
+ use mmlb4_m
+ implicit none
+ contains
+
+ subroutine MlbwCrossCalc_calculateCross(calculator, ener, Ipoten)
+ class(MlbwCrossCalc) :: calculator
+ real(kind=8)::ener
+ integer::ipoten
+
+ real(kind=8)::val
+ integer::i, itot
+
+ itot = calculator%covariance%getNumTotalParam()
+ call calculator%crossData%reserveColumns(calculator%row, itot + 1)
+ if (ener.eq.0.0d0) return
+
+ calculator%ener = ener
+ IF (calculator%ener.LT.0.0d0) calculator%ener = -calculator%ener
+ calculator%enerSq = dSQRT(calculator%ener)
+
+ CALL Abpart_Mlb (calculator)
+ CALL Parsh_Mlb ( calculator, Ipoten)
+
+
+ if (ener.lt.0.0d0) then
+ do i = 0, itot
+ val = calculator%crossData%getData(calculator%row, i, 1)
+ if (val.ne.0.0d0) then
+ ! addData accumulates. So to change the sign,
+ ! subtract twice the original value
+ val = -2.0d0*val
+ call calculator%crossData%addData(calculator%row, i, 1, val)
+ end if
+ end do
+ end if
+ end subroutine
+end module mmlb2_m
diff --git a/sammy/src/mlb/mmlb3.f b/sammy/src/mlb/mmlb3.f
deleted file mode 100644
index 632c5737c..000000000
--- a/sammy/src/mlb/mmlb3.f
+++ /dev/null
@@ -1,227 +0,0 @@
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Elastc_Mlb (Goj, Nent, Next, Zke,
- * If_Zkfe, Aaaone, Aaatwo, Aaathr,
- * Cs, Si, Dphi, Sigxxx, Dasigx, Iso, igr)
-C
-C *** Purpose -- Generate Sigxxx = elastic cross section
-C *** and Dasigx = partial derivative of elastic cross
-C *** section with respect to U-parameter
-C
- use fixedi_m
- use ifwrit_m
- use varyr_common_m
- use constn_common_m
- use templc_common_m, only : A_Ibr, A_Ibi, A_Ipr
- use EndfData_common_m, only : resparData, radFitFlags
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
- DIMENSION Zke(*), Aaaone(*), Aaatwo(*),
- * Aaathr(*),
- * Cs(*), Si(*), Dphi(*), Sigxxx(*), Dasigx(Nnnsig,*)
- DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
-C
-C
- Aaatot = Zero
- A = Zero
- X = Zero
- Sum = Zero
- Ntot = Nent + Next
- C = Goj*Pi100/Squ
- Jj = 0
- DO J=1,Nent
- Jj = Jj + J
- Aaatot = Aaatwo(J)
- DO I=1,Ntot
- IF (I.LE.J) THEN
- Ij = (J*(J-1))/2 + I
- ELSE
- Ij = (I*(I-1))/2 + J
- END IF
- Aaatot = Aaathr(Ij) + Aaatot
- END DO
- A = Two - Aaatot
- B = Two*Aaaone(J)
- E = Aaaone(J)**2 + 0.25D0*Aaatot**2
- D = One/Zke(J)**2
- IF (Krmatx.EQ.1) Sum = Sum + ( (One-Cs(J))*A + Si(J)*B + E ) *D
- IF (Krmatx.EQ.-1) Sum = Sum + ( (One-Cs(J))*A + Si(J)*B +
- * Aaathr(JJ) ) * D
- Ktru = radFitFlags%getTrueFitFlag(Igr, J)
- IF (Ktrue.gt.0.and.Ktru.LE.Ndasig) THEN
-cq Dasigx(1,Ktru) = Dasigx(1,Ktru) + X
-Cq *** note this is not right! don't have x defined yet! someday do it right!
- END IF
- Keff = radFitFlags%getEffFitFlag(Igr, J)
- IF (Keff.Gt.0.and.Keff.LE.Ndasig) THEN
- Dasigx(1,Keff) = Dasigx(1,Keff) +
- * ( (Si(J)*A+Cs(J)*B)*Dphi(J) ) * C / Zke(J)
- END IF
- END DO
- C = Goj*Pi100/Su
- Sigxxx(1) = C*Sum + Sigxxx(1)
- IF (Ksolve.EQ.2) RETURN
- Kiso = If_Zke
- IF (Kiso.gt.0.and.Kiso.le.Ndasig) THEN
- Dasigx(1,Kiso) = Dasigx(1,Kiso) + A/VarAbn
- END IF
-C
-C
- IF (Npr.GT.0) THEN
- DO Mm=1,Npr
- M = Kstart + Mm
- Sum = Zero
- Jj = 0
- DO J=1,Nent
- Jj = Jj + J
- Paatot = A_Ibi(J,M)
- DO I=1,Ntot
- IF (I.LE.J) THEN
- Ij = (J*(J-1))/2 + I
- ELSE
- Ij = (I*(I-1))/2 + J
- END IF
- Paatot = Paatot + A_Ipr(Ij,M)
- END DO
- IF (Krmatx.EQ.1) THEN
- A = -One + Cs(J) + 0.5*Aaatot
- B = Two*(SI(J)+Aaaone(J))
- Sum = Sum + (A*Paatot+B*A_Ibr(J,M))/Zke(J)**2
- ELSE IF (Krmatx.EQ.-1) THEN
- A = -One + Cs(J)
- B = Two*Si(J)
- Sum = Sum + (A*Paatot+B*A_Ibr(J,M)+A_Ipr(Jj,M)) /
- * Zke(J)**2
- ELSE
- STOP '[STOP in Elastc_Mlb in mlb/mmlb3.f]'
- END IF
- END DO
- Dasigx(1,M) = C*Sum + Dasigx(1,M)
- END DO
- END IF
- RETURN
- END
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Reactn_Mlb (Goj, Nent, Next, Zke, If_Zke,
- * Aaathr, Sigxxx, Dasigx, Iso, igr)
-C
-C *** Purpose -- Generate Sigxxx = reaction cross section
-C *** and Dasigx = partial derivative of reaction cross
-C *** section with respect to U-parameter
-C
- use fixedi_m
- use ifwrit_m
- use varyr_common_m
- use constn_common_m
- use templc_common_m, only : A_Ipr
- use EndfData_common_m, only : resparData, radFitFlags
- use SammySpinGroupInfo_M
- use SammyChannelInfo_M
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
- DIMENSION Zke(*), Aaathr(*),
- * Sigxxx(*), Dasigx(Nnnsig,*)
-C DIMENSION Aaathr(ntriag)
- DATA Zero /0.0d0/
-C
- IF (Next.EQ.0) RETURN
- Nent1 = Nent + 1
- NextN = Nent + Next
- Sum = Zero
- Sum1 = Zero
- Sum2 = Zero
- A = Goj*Pi100/Su
- C = Zero
- DO J=1,Nent
- B = A/Zke(J)**2
- Ktru = radFitFlags%getTrueFitFlag(Igr, J)
- IF (Ktru.GT.0) Dasigx(1,Ktru) = Dasigx(1,Ktru) + C
- DO I=Nent1,NextN
- IJ = (I*(I-1))/2 + J
- IF (I.EQ.Nent+1) Sum1 = Sum1 + Aaathr(IJ)*B
- IF (I.EQ.Nent+2) Sum2 = Sum2 + Aaathr(IJ)*B
- IF (I.GT.Nent+2) STOP '[STOP in Reactn_Mlb in mlb/mmlb3.f]'
- END DO
- END DO
- Sum = Sum1 + Sum2
- Sigxxx(1) = Sum + Sigxxx(1)
- Sig1 = Sum1 + Sig1
- Sig2 = Sum2 + Sig2
-C
- IF (Ksolve.EQ.2) RETURN
- Kiso = IF_Zke
- IF (Kiso.GT.0) THEN
- Dasigx(1,Kiso) = Dasigx(1,Kiso) + Sum/VarAbn
- END IF
- IF (Npr.GT.0) THEN
- DO Mm=1,Npr
- M = Kstart + Mm
- Sum = Zero
- DO J=1,Nent
- B = A/Zke(J)**2
- DO I=Nent1,Nextn
- Ij = (I*(I-1))/2 + J
- Sum = Sum + A_Ipr(Ij,M)*B
- END DO
- END DO
- Dasigx(1,M) = Sum + Dasigx(1,M)
- END DO
- END IF
- RETURN
- END
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Captur_Mlb (Goj, Nent, Zke, If_Zke, Aaatwo,
- * Sigxxx, Dasigx, Iso, igr)
-C
-C *** Purpose -- Generate Sigxxx = capture cross section
-C *** and Dasigx = partial derivative of capture cross
-C *** section with respect to U-parameter
-C
- use fixedi_m
- use ifwrit_m
- use varyr_common_m
- use constn_common_m
- use templc_common_m, only : A_Ibi
- use EndfData_common_m, only : resparData, radFitFlags
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
- DIMENSION Zke(*), Aaatwo(*),
- * Sigxxx(*), Dasigx(Nnnsig,*)
- Data Zero /0.0d0/
-C
- Sum = Zero
- B = Zero
- DO I=1,Nent
- Sum = Sum + Aaatwo(I)/Zke(I)**2
- Ktru = radFitFlags%getTrueFitFlag(Igr, I)
- IF (Ktru.GT.0) THEN
- Dasigx(1,Ktru) = Dasigx(1,Ktru) + B
- END IF
- END DO
- A = Goj*Pi100/Su
- Sigxxx(1) = A*Sum + Sigxxx(1)
- IF (Ksolve.EQ.2) RETURN
- Kiso = If_Zke
- IF (Kiso.GT.0) THEN
- Dasigx(1,Kiso) = Dasigx(1,Kiso) + A*Sum/VarAbn
- END IF
- IF (Npr.GT.0) THEN
- DO Mm=1,Npr
- M = Kstart + Mm
- Sum = Zero
- DO I=1,Nent
- Sum = Sum + A_Ibi(I,M)/Zke(I)**2
- END DO
- Dasigx(1,M) = A*Sum + Dasigx(1,M)
- END DO
- END IF
- RETURN
- END
diff --git a/sammy/src/mlb/mmlb3.f90 b/sammy/src/mlb/mmlb3.f90
new file mode 100644
index 000000000..876c996c3
--- /dev/null
+++ b/sammy/src/mlb/mmlb3.f90
@@ -0,0 +1,242 @@
+module mmlb3_m
+use MlbwCrossCalc_M
+IMPLICIT NONE
+contains
+!
+!
+! --------------------------------------------------------------
+!
+ SUBROUTINE Elastc_Mlb (calculator, Goj, Nent, Next, &
+ If_Zke, Iso, igr, ipar)
+!
+! *** Purpose -- Generate elastic cross section
+! *** and partial derivative of elastic cross
+! *** section with respect to U-parameter
+!
+ use constn_common_m, only : Pi100
+!
+ type(MlbwCrossCalc)::calculator
+ real(kind=8),intent(in)::Goj
+ integer, intent(in)::Nent, Next, If_Zke, Iso, igr, ipar
+
+ real(kind=8)::val, Zero, One, Two
+ real(kind=8)::A, Aaatot, B, C, E, D, Paatot, Sum, X
+ integer::I, Ij, J, Jj, M, Mm, Ntot, iflr
+ integer::Keff, Ktru, Kiso
+ DATA Zero /0.0d0/, One /1.0d0/, Two /2.0d0/
+!
+!
+ Aaatot = Zero
+ A = Zero
+ X = Zero
+ Sum = Zero
+ Ntot = Nent + Next
+ C = Goj*Pi100/calculator%enerSq
+ Jj = 0
+ DO J=1,Nent
+ Jj = Jj + J
+ Aaatot = calculator%Aaatwo(J, igr)
+ DO I=1,Ntot
+ IF (I.LE.J) THEN
+ Ij = (J*(J-1))/2 + I
+ ELSE
+ Ij = (I*(I-1))/2 + J
+ END IF
+ Aaatot = calculator%Aaathr(Ij, igr) + Aaatot
+ END DO
+ A = Two - Aaatot
+ B = Two*calculator%Aaaone(J, igr)
+ E = calculator%Aaaone(J, igr)**2 + 0.25D0*Aaatot**2
+ D = One/calculator%Zke(J, Igr)**2
+ IF (calculator%wantMlbw) then
+ Sum = Sum + ( (One-calculator%Cs(J))*A + calculator%Si(J)*B + E ) *D
+ else
+ Sum = Sum + ( (One-calculator%Cs(J))*A + calculator%Si(J)*B + &
+ calculator%Aaathr(JJ, igr) ) * D
+ end if
+ Ktru = calculator%radiusData%getTrueFitFlag(Igr, J)
+ IF (Ktru.gt.0) THEN
+!q val = calculator%crossData%getData(calculator%row, Ktru, 1) + X
+!q call calculator%crossData%addData(calculator%row, Ktru, 1, val)
+!q *** note this is not right! don't have x defined yet! someday do it right!
+ END IF
+ Keff = calculator%radiusData%getEffFitFlag(Igr, J)
+ IF (Keff.Gt.0) THEN
+ val = ( (calculator%Si(J)*A+calculator%Cs(J)*B)*calculator%Dphi(J) ) * C / calculator%Zke(J, Igr)
+ call calculator%crossData%addData(calculator%row, Keff, 1, val) ! crossData accumulates
+ END IF
+ END DO
+ C = Goj*Pi100/calculator%ener
+ call calculator%crossData%addData(calculator%row, 0, 1, C*Sum) ! crossData accumulates
+ IF (.not.calculator%wantDerivs) RETURN
+ Kiso = If_Zke
+ IF (Kiso.gt.0) THEN
+ val = A/calculator%getAbundance(Igr)
+ call calculator%crossData%addData(calculator%row, Kiso, 1, val) ! crossData accumulates
+ END IF
+!
+!
+ IF (calculator%inumSize.GT.0) THEN
+ DO Mm=1,calculator%inumSize
+ M = ipar + Mm ! position in calculator%Paathr, etc
+ iflr = calculator%inum(Mm,1) ! position in covariance
+ if (.not.calculator%covariance%contributes(Iflr)) cycle
+ Sum = Zero
+ Jj = 0
+ DO J=1,Nent
+ Jj = Jj + J
+ Paatot = calculator%Paatwo(J,M)
+ DO I=1,Ntot
+ IF (I.LE.J) THEN
+ Ij = (J*(J-1))/2 + I
+ ELSE
+ Ij = (I*(I-1))/2 + J
+ END IF
+ Paatot = Paatot + calculator%Paathr(Ij,M)
+ END DO
+ IF (calculator%wantMlbw) THEN
+ A = -One + calculator%Cs(J) + 0.5*Aaatot
+ B = Two*(calculator%SI(J)+calculator%Aaaone(J, igr))
+ Sum = Sum + (A*Paatot+B*calculator%Paaone(J,M))/calculator%Zke(J, Igr)**2
+ ELSE
+ A = -One + calculator%Cs(J)
+ B = Two*calculator%Si(J)
+ Sum = Sum + (A*Paatot+B*calculator%Paaone(J,M)+calculator%Paathr(Jj,M)) / &
+ calculator%Zke(J, Igr)**2
+ END IF
+ END DO
+ val = C*Sum
+ call calculator%crossData%addData(calculator%row, iflr, 1, val) ! crossData accumulates
+ END DO
+ END IF
+ RETURN
+ END
+!
+!
+! --------------------------------------------------------------
+!
+ SUBROUTINE Reactn_Mlb (calculator, Goj, Nent, Next, If_Zke, igr, ipar)
+!
+! *** Purpose -- Generate reaction cross section
+! *** and partial derivative of reaction cross
+! *** section with respect to U-parameter
+!
+ use constn_common_m, only : Pi100
+ use SammySpinGroupInfo_M
+ use SammyChannelInfo_M
+!
+ real(kind=8),intent(in)::Goj
+ integer,intent(in)::Nent, Next, If_Zke, igr, ipar
+ type(MlbwCrossCalc)::calculator
+
+ real(kind=8)::val, Zero, Sum, Sum1, Sum2
+ real(kind=8)::A,B,C
+ integer::Nent1,NextN, I, Ij, Mm, M, J, iflr
+ integer::Kiso, Ktru
+ DATA Zero /0.0d0/
+!
+ IF (Next.EQ.0) RETURN
+ Nent1 = Nent + 1
+ NextN = Nent + Next
+ Sum = Zero
+ Sum1 = Zero
+ Sum2 = Zero
+ A = Goj*Pi100/calculator%ener
+ C = Zero
+ DO J=1,Nent
+ B = A/calculator%Zke(J, Igr)**2
+ Ktru = calculator%radiusData%getTrueFitFlag(Igr, J)
+ IF (Ktru.GT.0) then
+ call calculator%crossData%addData(calculator%row, Ktru, 1, C)
+ end if
+ DO I=Nent1,NextN
+ IJ = (I*(I-1))/2 + J
+ IF (I.EQ.Nent+1) Sum1 = Sum1 + calculator%Aaathr(IJ, Igr)*B
+ IF (I.EQ.Nent+2) Sum2 = Sum2 + calculator%Aaathr(IJ, Igr)*B
+ IF (I.GT.Nent+2) STOP '[STOP in Reactn_Mlb in mlb/mmlb3.f]'
+ END DO
+ END DO
+ Sum = Sum1 + Sum2
+ call calculator%crossData%addData(calculator%row, 0, 1, Sum)
+!
+ IF (.not.calculator%wantDerivs) RETURN
+ Kiso = IF_Zke
+ IF (Kiso.GT.0) THEN
+ val = Sum/calculator%getAbundance(Igr)
+ call calculator%crossData%addData(calculator%row, Kiso, 1, val)
+ END IF
+ IF (calculator%inumSize.GT.0) THEN
+ DO Mm=1,calculator%inumSize
+ M = ipar + Mm ! position in calculator%Paathr, etc
+ iflr = calculator%inum(Mm,1) ! position in covariance
+ if (.not.calculator%covariance%contributes(Iflr)) cycle
+ Sum = Zero
+ DO J=1,Nent
+ B = A/calculator%Zke(J, Igr)**2
+ DO I=Nent1,Nextn
+ Ij = (I*(I-1))/2 + J
+ Sum = Sum + calculator%Paathr(Ij,M)*B
+ END DO
+ END DO
+ call calculator%crossData%addData(calculator%row, iflr, 1, Sum) ! crossData accumlates
+ END DO
+ END IF
+ RETURN
+ END
+!
+!
+! --------------------------------------------------------------
+!
+ SUBROUTINE Captur_Mlb (calculator, spinInfo, igr, ipar)
+!
+! *** Purpose -- Generate capture cross section
+! *** and partial derivative of capture cross
+! *** section with respect to U-parameter
+!
+ use constn_common_m, only : Pi100
+ use SammySpinGroupInfo_M
+!
+ type(MlbwCrossCalc)::calculator
+ type(SammySpinGroupInfo)::spinInfo
+ integer::imin, ipar
+
+ real(kind=8)::val, Zero, Sum, B, A, Goj
+ integer::I, M, Mm, Ktru, Kiso, igr, If_Zke, iflr
+ Data Zero /0.0d0/
+!
+ Sum = Zero
+ B = Zero
+ If_Zke = spinInfo%getAbundanceFitFlag()
+ Goj = calculator%getAbundance(igr)*spinInfo%getGFactor()
+ DO I=1,spinInfo%getNumEntryChannels()
+ Sum = Sum + calculator%Aaatwo(I, igr)/calculator%Zke(I, Igr)**2
+ Ktru = calculator%radiusData%getTrueFitFlag(Igr, I)
+ IF (Ktru.GT.0) THEN
+ call calculator%crossData%addData(calculator%row, Ktru, 1, B)
+ END IF
+ END DO
+ A = Goj*Pi100/calculator%ener
+ val = A*Sum
+ call calculator%crossData%addData(calculator%row, 0, 1, val)
+ IF (.not.calculator%wantDerivs) RETURN
+ Kiso = If_Zke
+ IF (Kiso.GT.0) THEN
+ val = A*Sum/calculator%getAbundance(Igr)
+ call calculator%crossData%addData(calculator%row, Kiso, 1, val)
+ END IF
+
+ IF (calculator%inumSize.GT.0) THEN
+ DO Mm=1,calculator%inumSize
+ M = ipar + Mm ! position in calculator%Paathr, etc
+ iflr = calculator%inum(Mm,1) ! position in covariance
+ if (.not.calculator%covariance%contributes(Iflr)) cycle
+ Sum = Zero
+ DO I=1,spinInfo%getNumEntryChannels()
+ Sum = Sum + calculator%Paatwo(I,M)/calculator%Zke(I, Igr)**2
+ END DO
+ call calculator%crossData%addData(calculator%row, Iflr, 1, A*Sum) ! crossData accumlates
+ END DO
+ END IF
+ RETURN
+ END
+end module mmlb3_m
diff --git a/sammy/src/mlb/mmlb4.f b/sammy/src/mlb/mmlb4.f
deleted file mode 100644
index ba662e229..000000000
--- a/sammy/src/mlb/mmlb4.f
+++ /dev/null
@@ -1,352 +0,0 @@
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Abpart_Mlb (Echan, Zke, Aaaone,
- * Aaatwo, Aaathr, P, Gami , Dphi)
-C
-C *** Purpose -- Generate Aaaone(I ,ig) = Gamma(I)*(E-Pken)/Den
-C *** Aaatwo(I ,ig) = Gamma(I)*Gamgam /Den
-C *** Aaathr(Ij,ig) = Gamma(I)*Gamma(J)/Den
-C *** and Paa's which are derivatives thereof wrt Iipar
-C *** Note that each term is summed over all resonances in group
-C
- use fixedi_m
- use ifwrit_m
- use fixedr_m
- use broad_common_m
- use varyr_common_m
- use SammyRMatrixParameters_M
- use SammyResonanceInfo_M
- use SammySpinGroupInfo_M
- use EndfData_common_m
- use RMatResonanceParam_M
- use xxx3
- use templc_common_m, only : I_Inotu, A_Ibr, A_Ibi, A_Ipr
- IMPLICIT none
-C
- type(SammyResonanceInfo)::resInfo
- type(SammySpinGroupInfo)::spinInfo
- type(RMatResonance)::resonance
- type(SammyParticlePairInfo)::pairInfo
- type(RMatParticlePair)::pair
- type(RMatChannelParams)::channel
- type(SammyChannelInfo)::channelInfo
- real(kind=8):: Echan(Ntotc,*), Zke(Ntotc,*),
- * Aaaone(Ntotc,*), Aaatwo(Ntotc,*), Aaathr(Ntriag,*),
- * P(*), Gami(*), Dphi(*)
- real(kind=8)::Zero, One
- real(kind=8)::D, De, Den, E, eres, G, Gamgam, Gamj
- real(kind=8)::gammaRed, Gamtot
- integer::I, ichan, Igr, Ij, Ipar, Ires, istart, Iflr
- integer::J, K, L, Min, Ntott, Ig
- real(kind=8)::W, Rho, width, X, Y, Z
-C
-C DIMENSION Ntot(Ngroup), Echan(Ntotc,Ngroup),
-C * Zke(ntotc,Ngroup),
-C * Aaaone(Ntotc,Ngroup),
-C * Aaatwo(Ntotc,Ngroup), Aaathr(Ntriag,Ngroup),
-C * P(Ntotc), Gami(Ntotc), Dphi(ntotc)
-C
- Data Zero /0.0d0/, One /1.0d0/
-C
- A_Ibr = 0.0d0
- A_Ibi = 0.0d0
- A_Ipr = 0.0d0
- CALL Zero_Array (Aaaone, Ntotc*Ngroup)
- CALL Zero_Array (Aaatwo, Ntotc*Ngroup)
- CALL Zero_Array (Aaathr, Ntriag*Ngroup)
-C
- Ipar = 0
- istart = 1
-
- DO Ig=1,resParData%getNumSpinGroups()
- call resparData%getSpinGroupInfo(spinInfo, Ig)
- Nnnn = Ig
- Min = istart
- Ij = 0
- Ntott = spinInfo%getNumChannels()
- DO I=1,Ntott
- call spinInfo%getChannelInfo(channelInfo, I)
- call resParData%getParticlePairInfo(
- * pairInfo,
- * channelInfo%getParticlePairIndex())
- call resParData%getParticlePair(pair, pairInfo)
- call resParData%getChannel(channel, channelInfo)
-
- L = channel%getL()
- Rho = channel%getApt()*Zke(I,Ig)*
- * dSQRT(Su-Echan(I,Nnnn))
- P(I) = One
- Dphi(I) = Zero
- IF (pair%getPnt().EQ.1) P(I) = Pfd (Rho, Dphi(I), L)
- END DO
-C
- DO Ires = Min, resParData%getNumResonances()
- call resParData%getResonanceInfo(resInfo, ires)
- istart = ires
- igr = resInfo%getSpinGroupIndex()
- if( igr.ne.ig) exit
-
- IF (resInfo%getIncludeInCalc()) THEN
- call resParData%getRedResonance(resonance, resInfo)
- ichan = spinInfo%getGammaWidthIndex()
- gammaRed = resonance%getWidth(ichan)
- De = Su - resonance%getEres()
- G = gammaRed**2
- DO I=1,Ntott
- ichan = spinInfo%getWidthForChannel(I)
- width = resonance%getWidth(ichan)
- G = G + P(I)*width**2
- END DO
- Gamtot = 2.0D0*G
- Den = De**2 + G**2
-C
- Gamgam = 2.0D0*gammaRed**2
- Ij = 0
- DO I=1,Ntott
- ichan = spinInfo%getWidthForChannel(I)
- width = resonance%getWidth(ichan)
- Gami(I) = 2.0*P(I)*width**2
- Aaaone(I,Ig) = Aaaone(I,Ig) + Gami(I)*De/Den
- Aaatwo(I,Ig) = Aaatwo(I,Ig) + Gami(I)*Gamgam/Den
- DO J=1,I
- Ij = Ij + 1
- ichan = spinInfo%getWidthForChannel(J)
- width = resonance%getWidth(ichan)
- Gamj = 2.0D0*P(J)*width**2
- Aaathr(Ij,Ig) = Aaathr(Ij,Ig) + Gami(I)*Gamj/Den
- END DO
- END DO
-C
- IF (Ksolve.NE.2) THEN
-C *** Note that, in what follows, X is prtl( (E-Pken)/Den ) wrt U-variable,
-C *** Y Gamgam /Den
-C *** and Z Gamma(J)/Den
-C
- IF (resInfo%getIncludeInCalc().and.
- * spinInfo%getIncludeInCalc()) THEN
-C
- IF (resInfo%getEnergyFitOption().gt.0) THEN
- Iflr = resInfo%getEnergyFitOption()
- Ipar = Ipar + 1
-C *** HERE RESONANCE ENERGY IS A VARIABLE
- IF (covData%contributes(Iflr)) THEN
- if (Iflr.ne.abs(I_Inotu(Ipar)))then
- STOP 'Count of varied resonance mmlb4'
- end if
- eres = resonance%getEres()
- X = 2.0D0*dSQRT(dABS(eres))*
- * (-One+2.0D0*De**2/Den)/Den
- Z = 4.0D0*dSQRT(dABS(eres))*De/Den**2
- Y = Z*Gamgam
- Ij = 0
- DO I=1,Ntott
- A_Ibr(I,Ipar) = Gami(I)*X
- A_Ibi(I,Ipar) = Gami(I)*Y
- DO J=1,I
- Ij = Ij + 1
- A_Ipr(IJ,Ipar) = Gami(I)*Gami(J)*Z
- END DO
- END DO
- END IF
- END IF
-C
- IF (resInfo%getChannelFitOption(1).NE.0) THEN
-C *** HERE Gamma-SUB-Gamma IS A VARIABLE
- Iflr = resInfo%getChannelFitOption(1)
- Ipar = Ipar + 1
- IF (covData%contributes(Iflr)) THEN
- if (Iflr.ne.abs(I_Inotu(Ipar)))then
- STOP 'Count of varied resonance mmlb4'
- end if
- X = -2.0D0*gammaRed*(De/Den)*
- * (Gamtot/Den)
- Y = 2.0D0*gammaRed*
- * (2.0-Gamgam*Gamtot/Den)/ Den
- Z = -2.0D0*gammaRed*Gamtot/Den**2
- Ij = 0
- DO I=1,Ntott
- A_Ibr(I,Ipar) = X*Gami(I)
- A_Ibi(I,Ipar) = Y*Gami(I)
- DO J=1,I
- Ij = Ij + 1
- A_Ipr(Ij,Ipar) =Z*Gami(I)*Gami(J)
- END DO
- END DO
- END IF
- END IF
-C
- DO K=1,Ntott
- IF (resInfo%getChannelFitOption(K+1).NE.0) THEN
-C
-C *** HERE Gamma-CHANNEL-K IS A VARIABLE
- Iflr = resInfo%getChannelFitOption(K+1)
- Ipar = Ipar + 1
- IF (covData%contributes(Iflr)) THEN
- if (Iflr.ne.abs(I_Inotu(Ipar)))then
- STOP 'Count of varied resonance mmlb4'
- end if
- ichan = spinInfo%getWidthForChannel(K)
- width = resonance%getWidth(ichan)
- W = 0.5D0*Gamtot*Gami(K)/Den
- X = 2.0D0*De/(width*Den)
- Y = 2.0D0*Gamgam/(width*Den)
- Z = 2.0D0/(width*Den)
- Ij = 0
- DO I=1,Ntott
- D = Zero
- IF (I.EQ.K) D = One
- A_Ibr(I,Ipar) = X*Gami(I)*(D-W)
- A_Ibi(I,Ipar) = Y*Gami(I)*(D-W)
- DO J=1,I
- Ij = Ij + 1
- E = Zero
- IF (J.EQ.K) E = One
- A_Ipr(Ij,Ipar) =
- * Z*Gami(I)*Gami(J)*(D+E-W)
- END DO
- END DO
- END IF
- END IF
- END DO
-C
- END IF
- END IF
- END IF
- END DO
- END DO
- RETURN
- END
-C
-C
-C --------------------------------------------------------------
-C
- SUBROUTINE Parsh_Mlb (
- * Zke,
- * Aaaone, Aaatwo, Aaathr, Pieces, Ipoten, Sigxxx, Dasigx,
- * Dbsigx)
-C
-C *** Purpose -- Form the cross section Sigxxx and the ( partial
-C *** derivatives of the cross section with respect to
-C *** the varied parameters ) = Dasigx
-C
- use over_common_m
- use oops_common_m
- use fixedi_m
- use ifwrit_m
- use exploc_common_m
- use fixedr_m
- use broad_common_m
- use varyr_common_m
- use EndfData_common_m
- use SammySpinGroupInfo_M
- use templc_common_m
- use xxx6
- IMPLICIT DOUBLE PRECISION (a-h,o-z)
-C
-C
- DIMENSION
- * Zke(Ntotc,*),
- * Aaaone(Ntotc,*),
- * Aaatwo(Ntotc,*), Aaathr(Ntriag,*), Pieces(*), Sigxxx(*),
- * Dasigx(*), Dbsigx(*)
-
- type(SammySpinGroupInfo)::spinInfo
- type(SammyResonanceInfo)::resInfo
-C
-C DIMENSION
-C * Aaaone(Ntotc,Ngroup), Aaatwo(Ntotc,Ngroup),
-C * Aaathr(ntriag,Ngroup), Cayiso(Numiso), Pilmis(Numiso),
-C * Pilsis(Numiso), Pieces(Ngroup)
-C
-C
-C
- NNF1 = 0
- NN2 = 0
- Kstart = 0
- Jstart = Nfpres
- imax = 0
-C
-C *** Do loop over groups (i.e., spin-parity groups) --
-C *** Goes to end of sub_routine
-C
- DO N=1,resParData%getNumSpinGroups()
- Kiso = 0
- call resParData%getSpinGroupInfo(spinInfo, N)
-
- imin = imax + 1
- Npr = 0
- DO I=imin,resParData%getNumResonances()
- call resParData%getResonanceInfo(resInfo, I)
- if( resInfo%getSpinGroupIndex().ne.n) exit
- imax = i
- IF (.not.resInfo%getIncludeInCalc()) cycle
-
-
- Mmax2 = spinInfo%getNumResPar()
- DO M=1,Mmax2
- if( m.eq.1) then
- Ifl = resInfo%getEnergyFitOption()
- else
- ifl = resInfo%getChannelFitOption(M-1)
- end if
- IF (Ifl.GT.0) THEN
- Npr = Npr + 1
- END IF
- end do
- END DO
-
- IF (spinInfo%getIncludeInCalc()) THEN
-C
- IF (Numiso.GT.0) THEN
- isoN = spinInfo%getIsotopeIndex()
- VarAbn = resParData%getAbundanceByIsotope(IsoN)
- Iso = isoN
- ELSE
- VarAbn = spinInfo%getAbundance()
- Iso = 1
- END IF
-C
- Nnnn = N
- Nnf1 = Nnf1 + Nn2
- Ntot = spinInfo%getNumChannels()
- Nn2 = Ntot*(Ntot+1)
- Nn = Nn2/2
-C
-C
-C *** CALCULATE SIN AND COS OF POTENTIAL SCATTERING PHASE
-C *** SHIFT, AND DERIVATIVE OF PHI WRT Rho
- IF (Kcros.LE.2) CALL Cossin (Zke(1,N),
- * A_Ics, A_Isi, A_Idphi, Nnnn, Ipoten)
-C
-C
- Agoj = varAbn*spinInfo%getGFactor()
-C *** SCATTERING CROSS SECTION (IE ELASTIC CROSS SECTION)
- nent = spinInfo%getNumEntryChannels()
- next = spinInfo%getNumExitChannels()
- iflAbund = spinInfo%getAbundanceFitFlag()
- IF (Kcros.EQ.1 .OR. Kcros.EQ.2) CALL Elastc_Mlb (Agoj,
- * Nent, Next, Zke(1,N), iflAbund,
- * Aaaone(1,N), Aaatwo(1,N), Aaathr(1,N),
- * A_Ics, A_Isi,
- * A_Idphi , Sigxxx, Dasigx, Iso, N)
-C
-C *** REACTION CROSS SECTIONS
- IF (Kcros.EQ.1 .OR. Kcros.EQ.3 .OR. Kcros.EQ.5) CALL
- * Reactn_Mlb (Agoj, Nent, Next, Zke(1,N), iflAbund,
- * Aaathr(1,N),
- * Sigxxx, Dasigx, Iso, N)
-C
-C *** CAPTURE CROSS SECTION
- IF (Kcros.EQ.1 .OR. Kcros.EQ.4 .OR. Kcros.EQ.5) CALL
- * Captur_Mlb (Agoj, Nent, Zke(1,N), iflAbund,
- * Aaatwo(1,N),
- * Sigxxx, Dasigx, Iso, N)
-C
- IF (Ksolve.NE.2) Kstart = Kstart + Npr
- IF (Kpiece.EQ.1) Pieces(N) = Sigxxx(1)
- END IF
- END DO
- RETURN
- END
diff --git a/sammy/src/mlb/mmlb4.f90 b/sammy/src/mlb/mmlb4.f90
new file mode 100644
index 000000000..20b99de51
--- /dev/null
+++ b/sammy/src/mlb/mmlb4.f90
@@ -0,0 +1,287 @@
+module mmlb4_m
+ use mmlb3_m
+ use MlbwCrossCalc_M
+ IMPLICIT none
+ contains
+!
+!
+! --------------------------------------------------------------
+!
+ SUBROUTINE Abpart_Mlb (calculator)
+!
+! *** Purpose -- Generate Aaaone(I ,ig) = Gamma(I)*(E-Pken)/Den
+! *** Aaatwo(I ,ig) = Gamma(I)*Gamgam /Den
+! *** Aaathr(Ij,ig) = Gamma(I)*Gamma(J)/Den
+! *** and Paa's which are derivatives thereof wrt Iipar
+! *** Note that each term is summed over all resonances in group
+!
+ use SammyRMatrixParameters_M
+ use SammyResonanceInfo_M
+ use SammySpinGroupInfo_M
+ use RMatResonanceParam_M
+ use xxx3
+!
+ type(MlbwCrossCalc)::calculator
+ type(SammyResonanceInfo)::resInfo
+ type(SammySpinGroupInfo)::spinInfo
+ type(RMatResonance)::resonance
+ type(SammyParticlePairInfo)::pairInfo
+ type(RMatParticlePair)::pair
+ type(RMatChannelParams)::channel
+ type(SammyChannelInfo)::channelInfo
+ real(kind=8)::Zero, One
+ real(kind=8)::D, De, Den, E, eres, G, Gamgam, Gamj
+ real(kind=8)::gammaRed, Gamtot
+ integer::I, ichan, Igr, Ij, Ipar, Ires, istart, Iflr, ii, m, ipos
+ integer::J, K, L, Min, Ntott, Ig
+ real(kind=8)::W, Rho, width, X, Y, Z
+
+ Data Zero /0.0d0/, One /1.0d0/
+!
+ IF (calculator%wantDerivs) then
+ calculator%Paaone = 0.0d0
+ calculator%Paatwo = 0.0d0
+ calculator%Paathr = 0.0d0
+ end if
+ calculator%Aaaone = 0.0d0
+ calculator%Aaatwo = 0.0d0
+ calculator%Aaathr = 0.0d0
+!
+ Ipar = 0
+ istart = 0
+
+ DO Ig=1,calculator%resData%getNumSpinGroups()
+ Min = istart + 1
+ call calculator%resData%getSpinGroupInfo(spinInfo, Ig)
+ call calculator%getParamPerSpinGroup(istart, ig)
+ Ij = 0
+ Ntott = spinInfo%getNumChannels()
+ DO I=1,Ntott
+ call spinInfo%getChannelInfo(channelInfo, I)
+ call calculator%resData%getParticlePairInfo( &
+ pairInfo, &
+ channelInfo%getParticlePairIndex())
+ call calculator%resData%getParticlePair(pair, pairInfo)
+ call calculator%resData%getChannel(channel, channelInfo)
+
+ L = channel%getL()
+
+ Rho = channel%getApt()*calculator%Zke(I,Ig)* &
+ dSQRT(calculator%ener-calculator%Echan(I,Ig))
+ calculator%Cs(I) = One
+ calculator%Dphi(I) = Zero
+ IF (pair%getPnt().EQ.1) calculator%Cs(I) = Pfd (Rho, calculator%Dphi(I), L)
+ END DO
+!
+ DO Ires = Min, calculator%resData%getNumResonances()
+ call calculator%resData%getResonanceInfo(resInfo, ires)
+ igr = resInfo%getSpinGroupIndex()
+ if( igr.ne.ig) exit
+
+ IF (resInfo%getIncludeInCalc()) THEN
+ call calculator%resData%getRedResonance(resonance, resInfo)
+ ichan = spinInfo%getGammaWidthIndex()
+ gammaRed = resonance%getWidth(ichan)
+ De = calculator%ener - resonance%getEres()
+ G = gammaRed**2
+ DO I=1,Ntott
+ ichan = spinInfo%getWidthForChannel(I)
+ width = resonance%getWidth(ichan)
+ G = G + calculator%Cs(I)*width**2
+ END DO
+ Gamtot = 2.0D0*G
+ Den = De**2 + G**2
+!
+ Gamgam = 2.0D0*gammaRed**2
+ Ij = 0
+ DO I=1,Ntott
+ ichan = spinInfo%getWidthForChannel(I)
+ width = resonance%getWidth(ichan)
+ calculator%Si(I) = 2.0*calculator%Cs(I)*width**2
+ calculator%Aaaone(I,Ig) = calculator%Aaaone(I,Ig) + calculator%Si(I)*De/Den
+ calculator%Aaatwo(I,Ig) = calculator%Aaatwo(I,Ig) + calculator%Si(I)*Gamgam/Den
+ DO J=1,I
+ Ij = Ij + 1
+ ichan = spinInfo%getWidthForChannel(J)
+ width = resonance%getWidth(ichan)
+ Gamj = 2.0D0*calculator%Cs(J)*width**2
+ calculator%Aaathr(Ij,Ig) = calculator%Aaathr(Ij,Ig) + calculator%Si(I)*Gamj/Den
+ END DO
+ END DO
+!
+ do ii = 1, calculator%inumSize
+! *** Note that, in what follows, X is prtl( (E-Pken)/Den ) wrt U-variable,
+! *** Y Gamgam /Den
+! *** and Z Gamma(J)/Den
+ if (calculator%inum(ii,3).ne.ires) cycle
+ if (.not.calculator%covariance%contributes( calculator%inum(ii,1))) cycle
+
+ ipos = ipar + ii
+ m = calculator%inum(ii,2)
+!
+
+!
+ IF (m.eq.1) THEN ! HERE RESONANCE ENERGY IS A VARIABLE
+
+ eres = resonance%getEres()
+ X = 2.0D0*dSQRT(dABS(eres))* &
+ (-One+2.0D0*De**2/Den)/Den
+ Z = 4.0D0*dSQRT(dABS(eres))*De/Den**2
+ Y = Z*Gamgam
+ Ij = 0
+ DO I=1,Ntott
+ calculator%Paaone(I,Ipos) = calculator%Si(I)*X
+ calculator%Paatwo(I,Ipos) = calculator%Si(I)*Y
+ DO J=1,I
+ Ij = Ij + 1
+ calculator%Paathr(IJ,Ipos) = calculator%Si(I)*calculator%Si(J)*Z
+ END DO
+ END DO
+
+
+!
+ else IF (m.eq.2) THEN ! HERE Gamma-SUB-Gamma IS A VARIABLE
+
+ X = -2.0D0*gammaRed*(De/Den)* &
+ (Gamtot/Den)
+ Y = 2.0D0*gammaRed* &
+ (2.0-Gamgam*Gamtot/Den)/ Den
+ Z = -2.0D0*gammaRed*Gamtot/Den**2
+ Ij = 0
+ DO I=1,Ntott
+ calculator%Paaone(I,Ipos) = X*calculator%Si(I)
+ calculator%Paatwo(I,Ipos) = Y*calculator%Si(I)
+ DO J=1,I
+ Ij = Ij + 1
+ calculator%Paathr(Ij,Ipos) =Z*calculator%Si(I)*calculator%Si(J)
+ END DO
+ END DO
+
+ else ! HERE Gamma-CHANNEL-K IS A VARIABLE
+ k = m - 2
+
+ ichan = spinInfo%getWidthForChannel(K)
+ width = resonance%getWidth(ichan)
+ W = 0.5D0*Gamtot*calculator%Si(K)/Den
+ X = 2.0D0*De/(width*Den)
+ Y = 2.0D0*Gamgam/(width*Den)
+ Z = 2.0D0/(width*Den)
+ Ij = 0
+ DO I=1,Ntott
+ D = Zero
+ IF (I.EQ.K) D = One
+ calculator%Paaone(I,Ipos) = X*calculator%Si(I)*(D-W)
+ calculator%Paatwo(I,Ipos) = Y*calculator%Si(I)*(D-W)
+ DO J=1,I
+ Ij = Ij + 1
+ E = Zero
+ IF (J.EQ.K) E = One
+ calculator%Paathr(Ij,Ipos) = &
+ Z*calculator%Si(I)*calculator%Si(J)*(D+E-W)
+ END DO
+ END DO
+
+
+ end if
+!
+
+ END Do
+ END IF
+ END DO
+ ipar = ipar + calculator%inumSize
+ END DO
+ RETURN
+ END
+!
+!
+! --------------------------------------------------------------
+!
+ SUBROUTINE Parsh_Mlb (calculator, Ipoten)
+!
+! *** Purpose -- Form the cross section and the ( partial
+! *** derivatives of the cross section with respect to
+! *** the varied parameters )
+!
+ use SammySpinGroupInfo_M
+ use xxx6
+!
+!
+ type(MlbwCrossCalc)::calculator
+ integer::Ipoten
+
+ type(SammySpinGroupInfo)::spinInfo
+ type(SammyResonanceInfo)::resInfo
+ real(kind=8)::Agoj
+ integer::I, Ifl, iflAbund, imax, Iso, M, N
+ integer::Nent, Next, NN2, NNF1, Ntot, imin, ipar
+!
+ NNF1 = 0
+ NN2 = 0
+
+ imax = 0
+ ipar = 0
+!
+! *** Do loop over groups (i.e., spin-parity groups) --
+! *** Goes to end of sub_routine
+!
+ DO N=1,calculator%resData%getNumSpinGroups()
+ call calculator%resData%getSpinGroupInfo(spinInfo, N)
+ imin = imax + 1
+ call calculator%getParamPerSpinGroup(imax, N)
+
+ IF (spinInfo%getIncludeInCalc()) THEN
+ Iso = spinInfo%getIsotopeIndex()
+!
+ Nnf1 = Nnf1 + Nn2
+ Ntot = spinInfo%getNumChannels()
+ Nn2 = Ntot*(Ntot+1)
+!
+!
+! *** CALCULATE SIN AND COS OF POTENTIAL SCATTERING PHASE
+! *** SHIFT, AND DERIVATIVE OF PHI WRT Rho
+ IF (calculator%reactType.LE.2) then
+ CALL Cossin (calculator%resData, &
+ calculator%Zke(1,N), &
+ calculator%Cs, &
+ calculator%Si, &
+ calculator%Dphi, &
+ N, &
+ Ipoten, &
+ calculator%enerSq, &
+ calculator%ener)
+ end if
+!
+!
+ Agoj = calculator%getAbundance(N)*spinInfo%getGFactor()
+! *** SCATTERING CROSS SECTION (IE ELASTIC CROSS SECTION)
+ nent = spinInfo%getNumEntryChannels()
+ next = spinInfo%getNumExitChannels()
+ iflAbund = spinInfo%getAbundanceFitFlag()
+ IF (calculator%reactType.EQ.1 .OR. &
+ calculator%reactType.EQ.2) then
+ CALL Elastc_Mlb (calculator, Agoj, Nent, Next, iflAbund, Iso, N, ipar)
+ end if
+!
+! *** REACTION CROSS SECTIONS
+ IF (calculator%reactType.EQ.1 .OR. &
+ calculator%reactType.EQ.3 .OR. &
+ calculator%reactType.EQ.5) then
+ CALL Reactn_Mlb (calculator, Agoj, Nent, Next, iflAbund, N, ipar)
+ end if
+!
+! *** CAPTURE CROSS SECTION
+ IF (calculator%reactType.EQ.1 .OR. &
+ calculator%reactType.EQ.4 .OR. &
+ calculator%reactType.EQ.5) then
+ CALL Captur_Mlb (calculator, spinInfo, N, ipar)
+ end if
+!
+ IF (calculator%wantCrossPerSpin) then
+ calculator%crossPerSpin(N) = calculator%crossData%getData(calculator%row, 0, 1)
+ end if
+ END IF
+ ipar = ipar + calculator%inumSize
+ END DO
+ RETURN
+ END
+end module mmlb4_m
diff --git a/sammy/src/sammy/CMakeLists.txt b/sammy/src/sammy/CMakeLists.txt
index bea4e3abd..3f7406096 100644
--- a/sammy/src/sammy/CMakeLists.txt
+++ b/sammy/src/sammy/CMakeLists.txt
@@ -240,11 +240,10 @@ APPEND_SET(SAMMY_SOURCES
../mas/mmasa.f
../mas/mmas6a.f
- ../mlb/mmlb0.f
- ../mlb/mmlb1.f
- ../mlb/mmlb2.f
- ../mlb/mmlb3.f
- ../mlb/mmlb4.f
+ ../mlb/mmlb0.f90
+ ../mlb/mmlb2.f90
+ ../mlb/mmlb3.f90
+ ../mlb/mmlb4.f90
../mlb/MlbwCrossCalc_M.f90
../mso/mmso0.f
diff --git a/sammy/src/the/CrossSectionCalcDriver_M.f90 b/sammy/src/the/CrossSectionCalcDriver_M.f90
index 37470210b..290c0faf6 100644
--- a/sammy/src/the/CrossSectionCalcDriver_M.f90
+++ b/sammy/src/the/CrossSectionCalcDriver_M.f90
@@ -136,6 +136,7 @@ module CrossSectionCalcDriver_M
end subroutine
subroutine CrossSectionCalcDriver_calcCross(this, row, ener, Ipoten)
+ use mmlb2_m
class(CrossSectionCalcDriver)::this
real(kind=8)::ener
integer::ipoten, row
@@ -149,6 +150,7 @@ module CrossSectionCalcDriver_M
if (this%mlbwInit) then
this%mlbw%row = row
+ call MlbwCrossCalc_calculateCross(this%mlbw, ener, Ipoten)
else if (this%artificalInit) then
this%artifical%row = row
call this%artifical%calculate(ener)
diff --git a/sammy/src/the/DerivativeHandler.cpp b/sammy/src/the/DerivativeHandler.cpp
index cdc52c745..13d1c9592 100644
--- a/sammy/src/the/DerivativeHandler.cpp
+++ b/sammy/src/the/DerivativeHandler.cpp
@@ -4,8 +4,8 @@ namespace sammy{
void DerivativeHandler::setUpList(SammyRMatrixParameters & params, AdjustedRadiusData & radii, int num){
int numIso = params.getNumIso();
- if (numIso < num) numIso = num;
- if (num == 1) numIso = 1; // numIso is sometimes redefined - why oh why
+ if (numIso < num) numIso = num; // this can happen if we want to keep track of channels included in the calculation with unique Echan values
+ if (num == 1) numIso = 1; // indicates that we can already sum and don't need to keep track of isotopes
if (numIso <= 1) {
if (getGridNumber() == 0) addGrid();
return;
@@ -25,6 +25,7 @@ namespace sammy{
if (!spinInfo->getIncludeInCalc()) includeInCalc = false;
int iso = spinInfo->getIsotopeIndex();
+ if (iso >= numIso) iso = numIso - 1;
// this allows legacy SAMMY to work. This shouldn't have been flagged
if (!includeInCalc) iso = 0;
diff --git a/sammy/src/the/mthe0.f90 b/sammy/src/the/mthe0.f90
index f03914621..e70e9c145 100644
--- a/sammy/src/the/mthe0.f90
+++ b/sammy/src/the/mthe0.f90
@@ -25,6 +25,7 @@ module mthe0_M
calcData, calcDataSelf, &
calcDataInit, calcDataSelfInit
use CrossSectionCalcDriver_M
+ use MlbwCalc_m
use xct_m
use SammyFlowControl_M, only : Where_To_Next
use cbro_common_m, only : Segnam
diff --git a/sammy/src/xxx/mxxx6.f90 b/sammy/src/xxx/mxxx6.f90
index b277de37f..448dda23d 100644
--- a/sammy/src/xxx/mxxx6.f90
+++ b/sammy/src/xxx/mxxx6.f90
@@ -3,22 +3,23 @@ module xxx6
!
! --------------------------------------------------------------
!
- SUBROUTINE Cossin (Zke, Cs, Si, Dphi, N, Ipoten)
+ SUBROUTINE Cossin (resParData, Zke, Cs, Si, Dphi, N, Ipoten, Squ, Su)
!
! *** PURPOSE -- GENERATE COS AND SIN OF ( 2.* Phi(HARD SHELL) )
!
- use sammy_CoulombSelector_I
- use fixedi_m
- use ifwrit_m
- use varyr_common_m
- use SammyResonanceInfo_M
- use EndfData_common_m
+ use sammy_CoulombSelector_I
+ use SammyResonanceInfo_M
use RMatResonanceParam_M
+ use SammySpinGroupInfo_M
+ use SammyChannelInfo_M
+ use SammyRMatrixParameters_M
implicit none
real(8), intent(in) :: Zke
real(8), intent(out) :: Cs, Si, Dphi
integer, intent(in) :: N, Ipoten
+ real(8), intent(in) :: Squ, Su
+ type(SammyRMatrixParameters)::resParData
real(8) :: A, A2, A4, B, C, D, Db, Dp, Ds, E, F, Fourpi, P, S, X, Y
integer :: iela, k, L, Nent
--
GitLab