Added some capability for LCOMP=2 format

openmc/data/endf.py

@@ -72,6 +72,24 @@ def float_endf(s):
     return float(ENDF_FLOAT_RE.sub(r'\1e\2', s))
 
 
+def int_endf(s):
+    """Conver string to int. Used for INTG records where blank entries
+    indicate a 0.
+
+    Parameters
+    ----------
+    s : str
+        Integer or spaces
+
+    Returns
+    -------
+    integer
+        The number or 0
+    """
+    s = s.strip()
+    return int(s) if s else 0
+
+
 def get_text_record(file_obj):
     """Return data from a TEXT record in an ENDF-6 file.
 
@@ -250,6 +268,50 @@ def get_tab2_record(file_obj):
 
     return params, Tabulated2D(breakpoints, interpolation)
 
+def get_intg_record(file_obj):
+    """
+    Return data from an INTG record in an ENDF-6 file.
+
+    Parameters
+    ----------
+    file_obj : file-like object
+        ENDF-6 file to read from
+
+    Returns
+    -------
+    array
+        The correlation matrix described in the INTG record
+    """
+    # determine how many items are in list and NDIGIT
+    items = get_cont_record(file_obj)
+    NDIGIT = int(items[2])
+    NNN = int(items[3]) # Number of parameters
+    NM = int(items[4]) # Lines to read
+    NROW_RULES = {2: 18,3: 12,4: 11,5: 9,6: 8}
+    NROW = NROW_RULES[NDIGIT]
+
+    # read lines and build correlation matrix
+    corr = np.identity(NNN)
+    for i in range(NM):
+        line = file_obj.readline()
+        ii = int_endf(line[:5]) - 1 #-1 to account for 0 indexing
+        jj = int_endf(line[5:10]) - 1
+        factor = 10**NDIGIT
+        for j in range(NROW):
+            if jj+j >= ii:
+                break
+            element = int_endf(line[11+(NDIGIT+1)*j:11+(NDIGIT+1)*(j+1)])
+            if element > 0:
+                corr[ii,jj] = (element+0.5)/factor
+            elif element < 0:
+                corr[ii,jj] = (element-0.5)/factor
+
+    #Symmetrize the correlation matrix
+    corr = corr + corr.T - np.diag(corr.diagonal())
+    return corr
+
+
+
 def get_evaluations(filename):
     """Return a list of all evaluations within an ENDF file.
 

openmc/data/resonance_covariance.py

@@ -6,7 +6,7 @@ from numpy.polynomial import Polynomial
 import pandas as pd
 
 from .data import NEUTRON_MASS
-from .endf import get_head_record, get_cont_record, get_tab1_record, get_list_record
+from .endf import get_head_record, get_cont_record, get_tab1_record, get_list_record, get_intg_record
 import openmc.checkvalue as cv
 from .resonance import ResonanceRange
 
@@ -80,7 +80,7 @@ class ResonanceCovariance(object):
                 formalism = items[3]  # resonance formalism
 
                 # Throw error for unsupported formalisms
-                if formalism in [0,1,7]:
+                if formalism in [0,7]:
                     raise TypeError('LRF= ', formalism,
                                     'covariance not supported for this formalism')
 
@@ -213,7 +213,44 @@ class MultiLevelBreitWignerCovariance(ResonanceRange):
 
             return mlbw
 
-        elif LCOMP in [0,2]:
+        elif LCOMP == 2: #Compact format - Resonances and individual
+                         #uncertainties followed by compact correlations
+            items, values = get_list_record(file_obj)
+            num_res = items[5] 
+            energy = values[0::12]
+            spin = values[1::12]
+            gt = values[2::12]
+            gn = values[3::12]
+            gg = values[4::12]
+            gf = values[5::12]
+            par_unc = []
+            for i in range(num_res):
+                res_unc = values[i*12+6:i*12+12]
+                #Delete 0 values (not provided in evaluation)
+                res_unc = [x for x in res_unc if x != 0.0]
+                par_unc.extend(res_unc)
+
+            records = []
+            for i, E in enumerate(energy):
+                records.append([energy[i], spin[i], gt[i], gn[i],
+                                gg[i], gf[i]])
+
+            corr = get_intg_record(file_obj)
+            cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc))
+
+            # Create pandas DataFrame with resonacne data
+            columns = ['energy', 'J', 'totalWidth', 'neutronWidth',
+                       'captureWidth', 'fissionWidth']
+            parameters = pd.DataFrame.from_records(records, columns=columns)
+
+            # Create instance of ReichMooreCovariance
+            mlbw = cls(energy_min, energy_max)
+            mlbw.parameters = parameters
+            mlbw.covariance = cov
+
+            return mlbw
+
+        elif LCOMP == 0 :
             raise TypeError('LCOMP = ' + str(LCOMP) + ' not supported')
 
 class SingleLevelBreitWignerCovariance(MultiLevelBreitWignerCovariance):
@@ -388,7 +425,44 @@ class ReichMooreCovariance(ResonanceRange):
 
             return rmc
 
-        elif LCOMP in [0,2]:
+        elif LCOMP == 2: #Compact format - Resonances and individual
+                         #uncertainties followed by compact correlations
+            items, values = get_list_record(file_obj)
+            num_res = items[5] 
+            energy = values[0::12]
+            spin = values[1::12]
+            gn = values[2::12]
+            gfa = values[3::12]
+            gfb = values[4::12]
+            par_unc = []
+            for i in range(num_res):
+                res_unc = values[i*12+6:i*12+12]
+                #Delete 0 values (not provided in evaluation)
+                res_unc = [x for x in res_unc if x != 0.0]
+                par_unc.extend(res_unc)
+
+            records = []
+            for i, E in enumerate(energy):
+                records.append([energy[i], spin[i], gn[i],
+                                gfa[i], gfb[i]])
+
+            corr = get_intg_record(file_obj)
+            cov = np.diag(par_unc).dot(corr).dot(np.diag(par_unc))
+
+            # Create pandas DataFrame with resonacne data
+            columns = ['energy', 'J', 'neutronWidth', 
+                       'fissionWidthA', 'fissionWidthB']
+            parameters = pd.DataFrame.from_records(records, columns=columns)
+
+            # Create instance of ReichMooreCovariance
+            rmc = cls(energy_min, energy_max)
+            rmc.parameters = parameters
+            rmc.covariance = cov
+
+            return rmc
+
+
+        elif LCOMP == 0:
             raise TypeError('LCOMP = ' + str(LCOMP) + ' not supported')