Bug in cccc.py when reading ISOTXS format

[valeriaRaffuzzi]

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!

[welcome]

Hi, and welcome to PyNE! 👋 Thanks for opening your first issue. We recommend that you include information such as the version of PyNE you're working with (eg, develop branch or a specific version), the platform you are operating on, the expected behavior, and the actual behavior you are bringing our attention to. The more deatil you provide, the better others in this community will be able to help you.