Bug in cccc.py when reading ISOTXS format

Describe the Bug
The function read in the class Isotxs in cccc.py to read multi group cross sections from the NJOY ISOTXS format returns wrong scattering matrices (at least when using the NJOY scattering matrix compacting flag IFOPT =1; I haven't tried with IFOPT = 2). Matrices are built by the function _read_nuclide_scatter, and stored inside nuc.micros["scat", block, g, fromg, order]. When trying to read from nuc.micros["scat",` block, g, fromg, order], the scattering matrices are wrong; the fix is to swap the loops for order in range(lordn): and for j in range(jl, ju + 1): inside _read_nuclide_scatter. That worked for me.

To Reproduce
Given an NJOY output tape25, formatted into ISOTXS format with flag IFOPT = 1, one can read it by:

# read NJOY output
isoFile = cccc.Isotxs('tape25')
isoFile.read()

# read isotope data
nuc = isoFile.find_nuclide('H1')

# save useful isotope data
ngroup  = isoFile.fc['ngroup']
blocks  = isoFile.fc['nscmax']

# initialise cross sections
P0scatt    = np.zeros((ngroup,ngroup,blocks))
P1scatt    = np.zeros((ngroup,ngroup,blocks))

# RETRIEVE CROSS SECTIONS 
for block in range(blocks):
    if (nuc.libParams['ords'][block] != 0):
        for gout in range(ngroup):
            for gin in range(ngroup):
                try:
                    P0scatt[gin,gout,block] = nuc.micros['scat',block,gout,gin,0]
                    P1scatt[gin,gout,block] = nuc.micros['scat',block,gout,gin,1]
                except KeyError:
                    pass

If then you plot P0scatt[gin,gout,block], where block corresponds to a different scattering reaction (0: inelastic, 1: elastic, 2: n2n, 3: total scattering), the matrix will look unphysical.

Screenshots or Code Snippets

The cccc.py code (lines approx. 360-380) could become:

for j in range(jl, ju + 1):
    g = j - 1
    assert g >= 0, "loading negative group in ISOTXS."
    jup = nuc.libParams["jj"][g, block] - 1
    jdown = nuc.libParams["jband"][g, block] - nuc.libParams["jj"][g, block]
    fromgroups = list(range(j - jdown, j + jup + 1))
    fromgroups.reverse()
    for order in range(lordn):
        for k in fromgroups:
            fromg = k - 1
            nuc.micros["scat", block, g, fromg, order] = r.get_float()[0]

NOTE: my solution works when IFOPT = 1, but I haven't tried with the NJOY flag IFOPT = 2. It should be tried before changing things!

pyne / pyne

Bug in cccc.py when reading ISOTXS format #1507