diff --git a/examples/create_mg_xs.sh b/examples/create_mg_xs.sh
index 9d4f0a4..01b72bf 100755
--- a/examples/create_mg_xs.sh
+++ b/examples/create_mg_xs.sh
@@ -63,3 +63,11 @@ cd $srcdir
 cd $scriptdir
 # The result should be a list of .data files in $SCRATCH_BARN/xs/pdtxs
 
+# Step 11: Run NJOY to generate the ACE files (can be done after step 2)
+cd $SCRATCH_BARN/xs
+./RunAce.sh
+# Step 12: Copy ACE files and create xsdir
+cd ace/xdata
+./copyAce.sh
+cd $scriptdir
+# The result should be ACE files and xsdir (tells MCNP where the ACE files are) in $SCRATCH_BARN/xs/ace/xdata
diff --git a/src/Readgroupr.py b/src/Readgroupr.py
index 4db94a0..aca97d4 100755
--- a/src/Readgroupr.py
+++ b/src/Readgroupr.py
@@ -200,6 +200,15 @@ def populate_data_dict(data, mfmtsSet, numGroups, numLegMoments, thermList, numS
             rxn['numDNGs'] = 6
         elif mf == 6 and mt == 18:
             rxn['numLegMoments'] = 1
+            rxn['flux'] = np.zeros((actualNumTherm, numGroups))
+            # rxn['lowEnergySpectrum'] = np.zeros((actualNumTherm, numGroups))
+            # rxnDict['lowEnergySpectrum'] = xsDict['lowEnergySpectrum'][::-1].copy()
+            # rxnDict['lowEnergyProd'] = xsDict['lowEnergyProd'][::-1].copy()
+            # rxnDict['highEnergyMatrix'] = xsDict['highEnergyMatrix'][::-1,::-1].copy()
+            # rxnDict['highestHighEnergyGroup'] = numGroups - xsDict['lowestHighEnergyGroup'] - 1
+            # rxnDict['promptChi'] = xsDict['promptChi'][::-1].copy()
+            # rxnDict['promptProd'] = xsDict['promptProd'][::-1].copy()
+            rxn['FissionMatrix'] = np.zeros((actualNumTherm, numGroups, numGroups))
         elif mf == 6:
             #'List' is indexed by temperature
             # Cannot set numLegMoments because some transfer matrices have fewer
@@ -340,15 +349,19 @@ def merge_data_prompt_fission_matrix(data, xsDict):
     mf,mt = xsDict['mf'], xsDict['mt']
     rxnDict = data['rxn'][(mf,mt)]
     numGroups = data['numGroups']
+    thermIndex = 0
+    if rxnDict['numTherm'] != 1:
+        thermIndex = data['thermList'].index(xsDict['temperature'])
+
     #The data is copied, not aliased
-    rxnDict['flux'] = xsDict['flux'][::-1].copy()
-    rxnDict['lowEnergySpectrum'] = xsDict['lowEnergySpectrum'][::-1].copy()
-    rxnDict['lowEnergyProd'] = xsDict['lowEnergyProd'][::-1].copy()
-    rxnDict['highEnergyMatrix'] = xsDict['highEnergyMatrix'][::-1,::-1].copy()
-    rxnDict['highestHighEnergyGroup'] = numGroups - xsDict['lowestHighEnergyGroup'] - 1
-    rxnDict['promptChi'] = xsDict['promptChi'][::-1].copy()
-    rxnDict['promptProd'] = xsDict['promptProd'][::-1].copy()
-    rxnDict['FissionMatrix'] = xsDict['FissionMatrix'][::-1, ::-1].copy()
+    rxnDict['flux'][thermIndex, :] = xsDict['flux'][::-1].copy()
+    # rxnDict['lowEnergySpectrum'] = xsDict['lowEnergySpectrum'][::-1].copy()
+    # rxnDict['lowEnergyProd'] = xsDict['lowEnergyProd'][::-1].copy()
+    # rxnDict['highEnergyMatrix'] = xsDict['highEnergyMatrix'][::-1,::-1].copy()
+    # rxnDict['highestHighEnergyGroup'] = numGroups - xsDict['lowestHighEnergyGroup'] - 1
+    # rxnDict['promptChi'] = xsDict['promptChi'][::-1].copy()
+    # rxnDict['promptProd'] = xsDict['promptProd'][::-1].copy()
+    rxnDict['FissionMatrix'][thermIndex, :, :] = xsDict['FissionMatrix'][::-1, ::-1].copy()
     #lowEnergySpectrum is indexed by (groupTo)
     #highEnergyMatrix is indexed by (groupFrom, groupTo)
     #highEnergyMatrix applies to groups 0 through highestHighEnergyGroup (inclusive)
@@ -1345,10 +1358,10 @@ def lookup_num_sig0(maxNumSig0, mf, mt):
     '''Look up the number of background xs the reaction contains (either maxNumSig0 or 1).'''
     if mf == 3:
         # NJOY 99:
-        mtsAtMultSig0 = {1, 2, 18, 102, 301, 318}
+        # mtsAtMultSig0 = {1, 2, 18, 102, 301, 318}
         # NJOY 2012:
-        #mtsAtMultSig0 = {1, 2, 18, 102, 301, 318}
-        #mtsAtMultSig0 |= {i for i in range(51,91+1)}
+        mtsAtMultSig0 = {1, 2, 18, 102, 301, 318}
+        mtsAtMultSig0 |= {i for i in range(51,91+1)}
         if mt in mtsAtMultSig0:
             return maxNumSig0
         else:
@@ -1372,7 +1385,7 @@ def lookup_num_therm(maxNumTherm, mf, mt):
         else:
             return 1
     elif mf == 6:
-        mtsDependOnT = [2]
+        mtsDependOnT = [2, 18]
         if (mt in mtsDependOnT) or (mt in thermalMTs):
             return maxNumTherm
         else:
@@ -1768,80 +1781,82 @@ def condense_xs(xsDataIn, energyMesh, flux, verbosity):
         numThermRxn = rxn['numTherm']
         numSig0Rxn = rxn['numSig0']
         numLegMomentsRxn = rxn['numLegMoments']
-        #
-        # First, we condense low-energy spectrum by column to an intermediate form (e <- g').
-        lowEnergySpectrumIn = rxnIn['lowEnergySpectrum']
-        lowEnergySpectrum = np.apply_along_axis(condenseFunc, 0, lowEnergySpectrumIn)
-        norm = np.sum(lowEnergySpectrum)
-        if norm:
-            lowEnergySpectrum /= norm
-        rxn['lowEnergySpectrum'] = lowEnergySpectrum
-        #
-        # Second, we condense prompt fission spectrum by column to an intermediate form (e <- g').
-        promptChiIn = rxnIn['promptChi']
-        promptChi = np.apply_along_axis(condenseFunc, 0, promptChiIn)
-        norm = np.sum(promptChi)
-        if norm:
-            promptChi /= norm
-        rxn['promptChi'] = promptChi
-        #
-        # Third, we condense prompt fission nu*sig_f (nu) using flux as weighting
-        promptProdIn = rxnIn['promptProd']
-        iTherm, iSig0 = 0, 0  # The 0th temperature and sig0 index is used for the weighting flux
-        wgt = flux[iTherm, iSig0, :] * promptProdIn
-        norm = elementFlux[iTherm, iSig0, :]
-        rxn['promptProd'] = np.bincount(energyMesh, weights=wgt) / norm
-        #
+        # #
+        # # First, we condense low-energy spectrum by column to an intermediate form (e <- g').
+        # lowEnergySpectrumIn = rxnIn['lowEnergySpectrum']
+        # lowEnergySpectrum = np.apply_along_axis(condenseFunc, 0, lowEnergySpectrumIn)
+        # norm = np.sum(lowEnergySpectrum)
+        # if norm:
+        #     lowEnergySpectrum /= norm
+        # rxn['lowEnergySpectrum'] = lowEnergySpectrum
+        # #
+        # # Second, we condense prompt fission spectrum by column to an intermediate form (e <- g').
+        # promptChiIn = rxnIn['promptChi']
+        # promptChi = np.apply_along_axis(condenseFunc, 0, promptChiIn)
+        # norm = np.sum(promptChi)
+        # if norm:
+        #     promptChi /= norm
+        # rxn['promptChi'] = promptChi
+        # #
+        # # Third, we condense prompt fission nu*sig_f (nu) using flux as weighting
+        # promptProdIn = rxnIn['promptProd']
+        # iTherm, iSig0 = 0, 0  # The 0th temperature and sig0 index is used for the weighting flux
+        # wgt = flux[iTherm, iSig0, :] * promptProdIn
+        # norm = elementFlux[iTherm, iSig0, :]
+        # rxn['promptProd'] = np.bincount(energyMesh, weights=wgt) / norm
+        # #
         # Fourth, we condense the full FissionMatrix
-        FissionMatrix = rxnIn['FissionMatrix']
-          # condense prompt fission matrix by column to an intermediate form (e <- g').
-        FissionMatrixIntermediate = np.apply_along_axis(condenseFunc, 1, FissionMatrix)
-        iTherm, iSig0 = 0, 0  # The 0th temperature and sig0 index is used for the weighting flux
-        wgt = np.dot(np.diag(flux[iTherm, iSig0, :]), FissionMatrixIntermediate[:,:])
-        norm = elementFlux[iTherm, iSig0, :]
-          # final condense for e <- e' (to element, from element).
-        rxn['FissionMatrix'] = np.dot(np.diag(1./norm), np.apply_along_axis(condenseFunc, 0, wgt))
-        #
-        # Fifth, we determine which elements are high-energy elements, if any
-        lastHighGroup = rxnIn['highestHighEnergyGroup']
-        if lastHighGroup == -1:
-            # The high-energy portion does not exist.
-            rxn['highestHighEnergyGroup'] = rxnIn['highestHighEnergyGroup']
-            rxn['highEnergyMatrix'] = rxnIn['highEnergyMatrix']
-        else:
-            highEnergyGroups = np.arange(0, lastHighGroup + 1)
-            highEnergyElements = np.unique(energyMesh[highEnergyGroups])
-            # We assume that the elements are intelligently ordered st low elements have higher energies
-            lastHighElement = highEnergyElements[-1]
-            rxn['highestHighEnergyGroup'] = lastHighElement
-            #
-            # We condense fission matrix by column to an intermediate form (e <- g').
-            highEnergyMatIn = rxnIn['highEnergyMatrix']
-            highEnergyMatIntermediate = np.apply_along_axis(condenseFunc, 1, highEnergyMatIn)
-            #
-            # Finally, we condense the fission matrix by row to the final form (e <- e').
-            # The 0th temperature and sig0 index is used for the weighting flux
-            rxn['highEnergyMatrix'] = np.zeros((lastHighElement+1, numElements))
-            highEnergyMat = rxn['highEnergyMatrix']
-            for elem in range(lastHighElement+1):
-                mask = (energyMesh == elem)
-                groupsInElem = np.where(energyMesh == elem)[0]
-                lowEnergyGroupsInElem = groupsInElem[groupsInElem > lastHighGroup]
-                highEnergyGroupsInElem = groupsInElem[groupsInElem <= lastHighGroup]
-                #
-                wgtHigh = flux[0, 0, highEnergyGroupsInElem]
-                wgtLow = flux[0, 0, lowEnergyGroupsInElem]
-                # The fission matrix in an element with high- and low-energy groups is a
-                # flux-weighted average
-                xsLow = xsData['rxn'][(3,18)]['xs'][0, 0, elem] * xsData['rxn'][(3,452)]['xs'][0, 0, elem]
-                sumLow = 0.
-                if xsLow:
-                    sumLow = xsLow * np.sum(wgtLow)
-                # The intermediate high energy fission matrix is nuSigf_{g'->e}
-                highEnergyMat[elem, :] = (
-                    np.sum(highEnergyMatIntermediate[highEnergyGroupsInElem, :] *
-                        wgtHigh[:, np.newaxis], axis=0) +
-                    sumLow * lowEnergySpectrum ) / (np.sum(wgtLow) + np.sum(wgtHigh))
+        rxn['FissionMatrix'] = np.zeros((numThermRxn, numElements, numElements))
+        for iTherm in range(numThermRxn):
+            FissionMatrix = rxnIn['FissionMatrix'][iTherm, :, :]
+              # condense prompt fission matrix by column to an intermediate form (e <- g').
+            FissionMatrixIntermediate = np.apply_along_axis(condenseFunc, 1, FissionMatrix)
+            # iTherm, iSig0 = 0, 0  # The 0th temperature and sig0 index is used for the weighting flux
+            wgt = np.dot(np.diag(flux[iTherm, iSig0, :]), FissionMatrixIntermediate[:,:])
+            norm = elementFlux[iTherm, iSig0, :]
+              # final condense for e <- e' (to element, from element).
+            rxn['FissionMatrix'][iTherm, :, :] = np.dot(np.diag(1./norm), np.apply_along_axis(condenseFunc, 0, wgt))
+        # #
+        # # Fifth, we determine which elements are high-energy elements, if any
+        # lastHighGroup = rxnIn['highestHighEnergyGroup']
+        # if lastHighGroup == -1:
+        #     # The high-energy portion does not exist.
+        #     rxn['highestHighEnergyGroup'] = rxnIn['highestHighEnergyGroup']
+        #     rxn['highEnergyMatrix'] = rxnIn['highEnergyMatrix']
+        # else:
+        #     highEnergyGroups = np.arange(0, lastHighGroup + 1)
+        #     highEnergyElements = np.unique(energyMesh[highEnergyGroups])
+        #     # We assume that the elements are intelligently ordered st low elements have higher energies
+        #     lastHighElement = highEnergyElements[-1]
+        #     rxn['highestHighEnergyGroup'] = lastHighElement
+        #     #
+        #     # We condense fission matrix by column to an intermediate form (e <- g').
+        #     highEnergyMatIn = rxnIn['highEnergyMatrix']
+        #     highEnergyMatIntermediate = np.apply_along_axis(condenseFunc, 1, highEnergyMatIn)
+        #     #
+        #     # Finally, we condense the fission matrix by row to the final form (e <- e').
+        #     # The 0th temperature and sig0 index is used for the weighting flux
+        #     rxn['highEnergyMatrix'] = np.zeros((lastHighElement+1, numElements))
+        #     highEnergyMat = rxn['highEnergyMatrix']
+        #     for elem in range(lastHighElement+1):
+        #         mask = (energyMesh == elem)
+        #         groupsInElem = np.where(energyMesh == elem)[0]
+        #         lowEnergyGroupsInElem = groupsInElem[groupsInElem > lastHighGroup]
+        #         highEnergyGroupsInElem = groupsInElem[groupsInElem <= lastHighGroup]
+        #         #
+        #         wgtHigh = flux[0, 0, highEnergyGroupsInElem]
+        #         wgtLow = flux[0, 0, lowEnergyGroupsInElem]
+        #         # The fission matrix in an element with high- and low-energy groups is a
+        #         # flux-weighted average
+        #         xsLow = xsData['rxn'][(3,18)]['xs'][0, 0, elem] * xsData['rxn'][(3,452)]['xs'][0, 0, elem]
+        #         sumLow = 0.
+        #         if xsLow:
+        #             sumLow = xsLow * np.sum(wgtLow)
+        #         # The intermediate high energy fission matrix is nuSigf_{g'->e}
+        #         highEnergyMat[elem, :] = (
+        #             np.sum(highEnergyMatIntermediate[highEnergyGroupsInElem, :] *
+        #                 wgtHigh[:, np.newaxis], axis=0) +
+        #             sumLow * lowEnergySpectrum ) / (np.sum(wgtLow) + np.sum(wgtHigh))
 
     return xsData
 
@@ -2041,21 +2056,24 @@ def interpolate_T_sig0_xs(data, desiredT, desiredSig0Vec, outputDict, verbosity=
             # Take a flux-weighted sum over all groups
             # Assumes you have (MF,MT) (3,18) (n,fission) and (3,452) (total nubar)
             flux = data['fluxOut']
-
-            fissionXS = data['rxn'][(3,18)]['xsOut']
-            nuBar = data['rxn'][(3,452)]['xsOut']
-            cutoffGroup = rxn['highestHighEnergyGroup']
-            lowEnergySpectrum = rxn['lowEnergySpectrum']
-            highEnergyMatrix = rxn['highEnergyMatrix']
-            effectiveSpectrum = np.zeros(numGroups)
-            effectiveSpectrum += lowEnergySpectrum * np.sum(
-                flux[cutoffGroup+1:] * fissionXS[cutoffGroup+1:] * nuBar[cutoffGroup+1:])
-            for groupFrom in range(0,cutoffGroup+1):
-                effectiveSpectrum += highEnergyMatrix[groupFrom, :] * flux[groupFrom]
-            effectiveSpectrum /= np.sum(effectiveSpectrum)
-            rxn['xsOut'] = effectiveSpectrum
-
-            FissionMatrix = rxn['FissionMatrix']
+            FissionMatrix = np.zeros((numGroups, numGroups))
+
+            # fissionXS = data['rxn'][(3,18)]['xsOut']
+            # nuBar = data['rxn'][(3,452)]['xsOut']
+            # cutoffGroup = rxn['highestHighEnergyGroup']
+            # lowEnergySpectrum = rxn['lowEnergySpectrum']
+            # highEnergyMatrix = rxn['highEnergyMatrix']
+            # effectiveSpectrum = np.zeros(numGroups)
+            # effectiveSpectrum += lowEnergySpectrum * np.sum(
+            #     flux[cutoffGroup+1:] * fissionXS[cutoffGroup+1:] * nuBar[cutoffGroup+1:])
+            # for groupFrom in range(0,cutoffGroup+1):
+            #     effectiveSpectrum += highEnergyMatrix[groupFrom, :] * flux[groupFrom]
+            # effectiveSpectrum /= np.sum(effectiveSpectrum)
+            # rxn['xsOut'] = effectiveSpectrum
+            for thermFrac, thermIndex in zip(thermFractions, thermIndices):
+                FissionMatrix += thermFrac * rxn['FissionMatrix'][thermIndex, :, :]
+            
+            rxn['FissionMatrixOut'] = FissionMatrix
             # Compute total fission source = F'*phi
             F_phi = np.dot(flux, FissionMatrix)
             # prompt Chi is normalized total fission source
@@ -2094,14 +2112,15 @@ def interpolate_T_sig0_xs(data, desiredT, desiredSig0Vec, outputDict, verbosity=
 
         flux = data['fluxOut']
         fission_xs = data['rxn'][(3 ,18)]['xsOut']
-        fission_x_prompt = data['rxn'][(6, 18)]['FissionMatrix']
+        fission_x_prompt = data['rxn'][(6, 18)]['FissionMatrixOut']
         promptProd = data['rxn'][(6, 18)]['promptProd']
         nu_prompt = promptProd/fission_xs
         nu_delayed = data['rxn'][(3 ,455)]['xsOut']
         chis_delayed = data['rxn'][(5 ,455)]['delayedChi']
         chi_delayed = np.sum(chis_delayed, axis=0)
 
-        nu_ss = nu_prompt + nu_delayed
+#debug       nu_ss = nu_prompt + nu_delayed
+        nu_ss = (nu_prompt + nu_delayed) * data['rxn'][(3, 18)]['xsOut']
         n_per_gout = ( np.dot(flux, fission_x_prompt) + \
                    chi_delayed*np.sum(nu_delayed*fission_xs*flux) )
         chi_ss = n_per_gout/np.sum(n_per_gout)
@@ -2210,9 +2229,15 @@ def write_pdt_xs(filePath, data, temperature, format='csr', whichXS='all', fromF
         if (3,18) in mfmts:
             # add MT 18
             numXS += 1
+        if (3, 452) in mfmts:
+            # add MT 452 (total nu)
+            numXS += 1
         if (3, 455) in mfmts:
             # add MT 455 (delayed nu)
             numXS += 1
+        if (3, 456) in mfmts:
+            # add MT 456 (prompt nu)
+            numXS += 1
         if (5,455) in mfmts:
             # Add decayConst and delayedChi for delay fission neutrons
             # Both are regarded as 1D xs (MF 3)
@@ -2287,8 +2312,12 @@ def write_pdt_xs(filePath, data, temperature, format='csr', whichXS='all', fromF
         mtsForMF3 = []
         if (3,18) in mfmts:
             mtsForMF3 = [18]
+        if (3,452) in mfmts:
+            mtsForMF3.append(452)
         if (3,455) in mfmts:
             mtsForMF3.append(455)
+        if (3,456) in mfmts:
+            mtsForMF3.append(456)
         if whichXS.lower().strip() in ['all', 'everything']:
             mtsForMF3 = [mt for (mf,mt) in sorted(mfmts) if (mf == 3 and mt != 1)]
         elif whichXS.lower().strip() in 'invel':
@@ -2320,7 +2349,7 @@ def write_pdt_xs(filePath, data, temperature, format='csr', whichXS='all', fromF
             fid.write(multiline_string(xs, 20, 5, 12))
             # write total fission transfer matrix (chi dot nusigf)
             pdtMT = 2518
-            xs = data['rxn'][(6,18)]['FissionMatrix']
+            xs = data['rxn'][(6,18)]['FissionMatrixOut']
             fid.write('MT {0}, Moment {1}\n'.format(pdtMT, 0))
             for g in range(numGroups):
                 fid.write('{0}, first, last: '.format('  Sink'))
diff --git a/src/indicators.py b/src/indicators.py
index 14ac4b7..d9283cb 100755
--- a/src/indicators.py
+++ b/src/indicators.py
@@ -207,9 +207,9 @@ def do_all(inputDict):
             materials.append(mat.get_triga_grid_plate_material())
             materials.append(mat.get_triga_lead_material())
     elif materialOpt == 'ctrigafuel':
-        materials.append(mat.get_depleted_triga_fuel_material())
+        materials.append(mat.get_ctriga_fuel_material())
         if workOpt == 'wgt':
-            materials.append(mat.get_triga_clad_material())
+            materials.append(mat.get_ctriga_clad_material())
             materials.append(mat.get_triga_zirconium_material())
             materials.append(mat.get_triga_graphite_material())
             materials.append(mat.get_triga_borated_graphite_material())
@@ -242,6 +242,19 @@ def do_all(inputDict):
             materials.append(mat.get_triga_air_material())
             materials.append(mat.get_triga_grid_plate_material())
             materials.append(mat.get_triga_lead_material())
+    elif materialOpt == 'CASL':
+        materials.append(mat.get_CASL_fuel_p4_211_material())
+        materials.append(mat.get_CASL_fuel_p4_262_material())
+        materials.append(mat.get_CASL_fuel_p5_31_material())
+        if workOpt == 'wgt':
+            materials.append(mat.get_CASL_cladding_p5_material())
+            materials.append(mat.get_CASL_pyrex_p5_material())
+            materials.append(mat.get_CASL_AIC_p5_material())
+            materials.append(mat.get_CASL_B4C_p5_material())
+            materials.append(mat.get_CASL_moderator_p5_material())
+            materials.append(mat.get_CASL_gas_p5_material())
+            materials.append(mat.get_CASL_coreplates_material())
+            materials.append(mat.get_CASL_StainlessSteel_p5_material())
     elif materialOpt == 'deb':
         materials.append(mat.get_bruss_enriched_rod_fuel_material())
     else:
@@ -1303,7 +1316,7 @@ def define_input_parser():
     parser.add_argument('-G', '--groupopt', help="Base coarse group structure file head. If the flux work option is used, the output group structure is the base coarse group structure with the RRR overwritten by the hyperfine group structure. Some examples include 'c5g7', 'scale-44', and 'res-N' where N=1,..9. Always looks in default directory of makegroups.", default=defaults['groupopt'])
     parser.add_argument('-r', '--resolution', help='Resolution to use for the pointwise flux calculations', type=int, choices=range(11), default=defaults['resolution'])
     parser.add_argument('-E', '--energyspacing', help='The maximum energy jump in the final grid is equal to this factor multiplied by a jump normalization based on downscattering distance off a heavy nucleus.', type=float, default=defaults['energyspacing'])
-    parser.add_argument('-m', '--materialopt', help=" Unless 'manual' is used, specifies a set of materials to use. If 'manual' is used, give a space-separated list of material names in 'listmaterials'.", choices=['3','4','5','c5g7', 'graphite', 'iron', 'kpin', 'kenrichedpin', 'kcladpin', 'kpin2d', 'kenrichedpin2d', 'kmoxpin2d', 'kmoxenrichedpin2d', 'trigafuel', 'ctrigafuel', 'ctrigafuel_0', 'trigamore', 'deb', 'manual'], default=defaults['materialopt'])
+    parser.add_argument('-m', '--materialopt', help=" Unless 'manual' is used, specifies a set of materials to use. If 'manual' is used, give a space-separated list of material names in 'listmaterials'.", choices=['3','4','5','c5g7', 'graphite', 'iron', 'kpin', 'kenrichedpin', 'kcladpin', 'kpin2d', 'kenrichedpin2d', 'kmoxpin2d', 'kmoxenrichedpin2d', 'trigafuel', 'ctrigafuel', 'ctrigafuel_0', 'trigamore', 'CASL', 'deb', 'manual'], default=defaults['materialopt'])
     parser.add_argument('-i', '--indicatormaterials', dest='listmaterials', help="When manual 'materialopt' is used, specify the materials to use.", nargs='+', default=defaults['listmaterials'], choices=mat.get_materials_name2function_dict().keys())
     parser.add_argument('-t', '--temperaturedependence', help="If specified, use the temperature in the PENDF file that corresponds to the temperature of the material (as specified in 'materials_materials.py'). If this temperature does not exist, an error is thrown. If not specified, the first temperature in the PENDF file is used. This flag is added for 'materialopt' of '3'.", action='store_true', default=False)
     parser.add_argument('-M', '--meshname', help="The energy mesh on which the output fluxes are constructed. {r} is replaced by the resolution. If the 'flux' workopt is used, the mesh file is written to. Else, is it read from. If no file type is specified, '.txt' is used. If no directory is specified, the 'energyDat' directory from directories is used.", default=defaults['meshname'])
diff --git a/src/indicators_clustering.py b/src/indicators_clustering.py
index 85beb87..a082fab 100755
--- a/src/indicators_clustering.py
+++ b/src/indicators_clustering.py
@@ -220,6 +220,8 @@ def read_energy_and_indicators(dataDict):
         importanceDict = {'tdFUEL_0': 1}
     elif materialOpt == 'trigamore':
         importanceDict = {'tFUEL': 10, 'tCLAD': 2, 'tZIRC': 2, 'tIRRADIATIONTUBE': 1}
+    elif materialOpt == 'CASL':
+        importanceDict = {'CASLfuel_P4_211': 1, 'CASLfuel_P4_262': 1, 'CASLfuel_P5_31': 1}
     elif materialOpt == 'manual':
         if not importancesList:
             importancesList = [1]*len(materialsList)
@@ -976,7 +978,7 @@ def define_input_parser():
     parser.add_argument('-n', '--numcoarsegroups', help="The number of coarse groups to be used. If nonzero, overwrites the internal members of 'coarsebdrs'", type=int, default=defaults['numcoarsegroups'])
     parser.add_argument('-l', '--listnumelements', help='Number of elements to be put in each coarse boundary. Number of arguments should be one less than the number of coarse boundaries. Takes priority over "elements" if set', type=int, nargs='+', default=defaults['listnumelements'])
     parser.add_argument('-r', '--resolution', help='Resolution to use for the pointwise flux calculations', type=int, choices=range(11), default=defaults['resolution'])
-    parser.add_argument('-m', '--materialopt', help="Unless 'manual' is used, specifies a set of materials to use. If 'manual' is used, give a space-separated list of material names in 'listmaterials'.", choices=['4','5','c5g7', 'graphite', 'iron', 'kpin', 'kenrichedpin', 'kcladpin', 'kpin2d', 'kenrichedpin2d', 'kmoxpin2d', 'kmoxenrichedpin2d', 'trigafuel', 'ctrigafuel', 'ctrigafuel_0' 'trigamore', 'manual'], default=defaults['materialopt'])
+    parser.add_argument('-m', '--materialopt', help="Unless 'manual' is used, specifies a set of materials to use. If 'manual' is used, give a space-separated list of material names in 'listmaterials'.", choices=['4','5','c5g7', 'graphite', 'iron', 'kpin', 'kenrichedpin', 'kcladpin', 'kpin2d', 'kenrichedpin2d', 'kmoxpin2d', 'kmoxenrichedpin2d', 'trigafuel', 'ctrigafuel', 'ctrigafuel_0', 'trigamore', 'CASL', 'manual'], default=defaults['materialopt'])
     parser.add_argument('-i', '--indicatormaterials', dest='listmaterials', help="When manual 'materialopt' is used, specify the materials to use.", nargs='+', default=defaults['listmaterials'])
     parser.add_argument('-I', '--importances', dest='listimportances', help="When manual 'materialopt' is used, specify the weightings (importances) to use when clustering.", nargs='+', type=int, default=[])
     parser.add_argument('-p', '--plotopt', help='Which observations to plot', choices=['none', 'first', 'last', 'firstlast', 'half', 'all', 'sum'], default=defaults['plotopt'])
diff --git a/src/materials_bondarenko.py b/src/materials_bondarenko.py
index ab0c93b..40f31cd 100755
--- a/src/materials_bondarenko.py
+++ b/src/materials_bondarenko.py
@@ -256,9 +256,9 @@ def form_and_print_macroscopic_xs(dirr, ZAList, material, numGroups, verbosity=F
     for MT in MTsAvg:
         if np.all(norms[MT]>0.0):
             xsOut[MT] /= norms[MT]
-
+    
     # Recompute steady-state nu and chi
-    if all(MTs in xsIn for MTs in [MTnuSigF, MTfission, MTssNu, MTssChi]):
+    if all(mts in MTs for mts in [MTnuSigF, MTfission, MTssNu, MTssChi]):
         flux = xsOut[MTwgt]
         promptProd = xsOut[MTnuSigF]
         fission_xs = xsOut[MTfission]
@@ -268,7 +268,7 @@ def form_and_print_macroscopic_xs(dirr, ZAList, material, numGroups, verbosity=F
         fission_x_prompt = xsOut[MTfissionMatrix]
 
         nu_prompt = promptProd/fission_xs
-        nu_ss = nu_prompt + nu_delayed
+        nu_ss = (nu_prompt + nu_delayed) * fission_xs
         n_per_gout = ( np.dot(fission_x_prompt, flux) + \
                        chi_delayed*np.sum(nu_delayed*fission_xs*flux) )
         chi_ss = n_per_gout/np.sum(n_per_gout)
diff --git a/src/materials_materials.py b/src/materials_materials.py
index 86a1d94..731ed8c 100644
--- a/src/materials_materials.py
+++ b/src/materials_materials.py
@@ -32,11 +32,37 @@ def get_materials_name2function_dict():
         'certHh2o': get_cert_Hh2o,
         # For CASL use
         'CASLfuel': get_CASL_fuel_material,
+        'CASLfuel_900K': get_CASL_fuel_material_900K,
+        'CASLfuel_1200K': get_CASL_fuel_material_1200K,
         'CASLgas': get_CASL_gas_material,
         'CASLclad': get_CASL_cladding_material,
         'CASLmod': get_CASL_moderator_material,
+        'CASLmod1B': get_CASL_moderator1B_material,
+        'CASLmod2E': get_CASL_moderator2E_material,
         'CASLss': get_CASL_StainlessSteel_material,
         'CASLpyrex': get_CASL_pyrex_material,
+        'CASLb4c': get_CASL_B4C_material,
+        'CASLaic': get_CASL_AIC_material,
+        # CASL Problem 4 (Page.54-55)
+        'CASLfuel_P4_211': get_CASL_fuel_p4_211_material,
+        'CASLfuel_P4_262': get_CASL_fuel_p4_262_material,
+        'CASLmod_P4': get_CASL_moderator_p4_material,
+        'CASLtopnozzle': get_CASL_topnozzle_material,
+        'CASLbottomnozzle': get_CASL_bottomnozzle_material,
+        'CASLcoreplates': get_CASL_coreplates_material,
+        # CASL Problem 5 (Page.62 & 67)
+        'CASLfuel_P5_31': get_CASL_fuel_p5_31_material,
+        'CASLmod_P5': get_CASL_moderator_p5_material,
+        'CASLpyrex_P5': get_CASL_pyrex_p5_material,
+        'CASLb4c_P5': get_CASL_B4C_p5_material,
+        'CASLaic_P5': get_CASL_AIC_p5_material,
+        'CASLss_P5': get_CASL_StainlessSteel_p5_material,
+        'CASLgas_P5': get_CASL_gas_p5_material,
+        'CASLclad_P5': get_CASL_cladding_p5_material,
+
+
+        # For 410 design use
+        '410Design': get_410Design_material,
         # Simple examples
         'hpu': get_hpu_slurry_material,
         'hheu': get_hheu_slurry_material,
@@ -67,15 +93,18 @@ def get_materials_name2function_dict():
         'cFCHAMBER': get_c5g7_fission_chamber_material,
         'cCR': get_c5g7_control_rod_material,
         # TRIGA (BOL)
-        'tcFUEL': get_depleted_triga_fuel_material,
+        'tdFUEL': get_depleted_triga_fuel_material,
         'tdFUEL_0': get_depleted_triga_fuel_material_0,
         'tFUEL': get_triga_fuel_material,
+        'tcFUEL': get_ctriga_fuel_material,
         'tZIRC': get_triga_zirconium_material,
         'tCLAD': get_triga_clad_material,
+        'tcCLAD': get_ctriga_clad_material,
         'tMOD': get_triga_moderator_material,
         'tAIR': get_triga_air_material,
         'tGRIDPLATE': get_triga_grid_plate_material,
         'tGRAPHITE': get_triga_graphite_material,
+        'tWATERGRAPHITE': get_triga_watergraphite_material,
         'tBORATEDGRAPHITE': get_triga_borated_graphite_material,
         'tB4C': get_triga_b4c_material,
         'tLEAD': get_triga_lead_material,
@@ -125,6 +154,9 @@ def get_materials_name2function_dict():
         'mtTGRAPHITE_9': get_multi_temperature_triga_graphite_material_T9,
         'mtDTFUEL_0': get_multi_temperature_depleted_triga_fuel_material_T0,
         'mtDTFUEL_7': get_multi_temperature_depleted_triga_fuel_material_T7,
+        'mtTB4C_0': get_multi_temperature_triga_b4c_material_T0,
+        'mtTB4C_5': get_multi_temperature_triga_b4c_material_T5,
+        'mtTB4C_7': get_multi_temperature_triga_b4c_material_T7,
     }
 
 ###############################################################################
@@ -151,6 +183,52 @@ def get_CASL_fuel_material():
         elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
     return material
 
+def get_CASL_fuel_material_900K():
+    shortName = 'CASLfuel_900K'
+    longName = 'CASL fuel material at 900K'
+    massDensity = 10.257 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 900. #K
+    thermalOpt = 'free'
+    uAtomFractionsDict = {234:6.11864E-06, 235:7.18132E-04, 236:03.29861E-06, 238:2.21546E-02}
+    elemAtomFracDict = {'O':4.57642E-02, 'U':2.28821E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, uAtomFractionsDict, 'U')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_fuel_material_1200K():
+    shortName = 'CASLfuel_1200K'
+    longName = 'CASL fuel material at 1200K'
+    massDensity = 10.257 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 1200. #K
+    thermalOpt = 'free'
+    uAtomFractionsDict = {234:6.11864E-06, 235:7.18132E-04, 236:03.29861E-06, 238:2.21546E-02}
+    elemAtomFracDict = {'O':4.57642E-02, 'U':2.28821E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, uAtomFractionsDict, 'U')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
 def get_CASL_gas_material():
     shortName = 'CASLgas'
     longName = 'CASL gas material'
@@ -212,6 +290,29 @@ def get_CASL_moderator_material():
     longName = 'CASL moderator material'
     massDensity = 0.743 #g/cc
     fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'h2o'
+    bAtomFractionsDict = {10:1.07070E-05, 11:4.30971E-05}
+    elemAtomFracDict = {'H':4.96224E-02, 'O':2.48112E-02, 'B':5.38041E-05}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_moderator2E_material():
+    shortName = 'CASLmod2E'
+    longName = 'CASL moderator 2E material'
+    massDensity = 0.743 #g/cc
+    fuelRadius = 0.4096 #cm
     temperature = 600. #K
     thermalOpt = 'h2o'
     bAtomFractionsDict = {10:1.07070E-05, 11:4.30971E-05}
@@ -230,6 +331,29 @@ def get_CASL_moderator_material():
         elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
     return material
 
+def get_CASL_moderator1B_material():
+    shortName = 'CASLmod1B'
+    longName = 'CASL moderator 1B material'
+    massDensity = 0.661 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 600. #K
+    thermalOpt = 'h2o'
+    bAtomFractionsDict = {10:9.52537E-06, 11:3.83408E-05}
+    elemAtomFracDict = {'H':4.41459E-02, 'O':2.20729E-02, 'B':4.786617E-05}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
 def get_CASL_StainlessSteel_material():
     shortName = 'CASLss'
     longName = 'CASL stainless-steel material'
@@ -285,6 +409,465 @@ def get_CASL_pyrex_material():
         elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
     return material
 
+def get_CASL_B4C_material():
+    shortName = 'CASLb4c'
+    longName = 'CASL B4C material'
+    atomDensity = 9.59100E-02 #1/barn-cm
+    fuelRadius = 0.4096 #cm
+    temperature = 600. #K
+    thermalOpt = 'free'
+    bAtomFractionsDict = {10:1.52689E-02, 11:6.14591E-02}
+    elemAtomFracDict = {'B':7.67280E-02, 'C':1.91820E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    override_abundances_as_elemental(ZAList, abundanceDict, 'C')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, atomDensity=atomDensity)
+    return material
+
+def get_CASL_AIC_material():
+    shortName = 'CASLaic'
+    longName = 'CASL AIC material'
+    atomDensity = 5.38440718E-02 #1/barn-cm
+    fuelRadius = 0.4096 #cm
+    temperature = 600. #K
+    thermalOpt = 'free'
+    agAtomFractionsDict = {107:2.36159E-02, 109:2.19403E-02}
+    cdAtomFractionsDict = {106:3.41523E-05, 108:2.43165E-05, 110:3.41250E-04, 111:3.49720E-04, 112:6.59276E-04, 113:3.33873E-04, 114:7.84957E-04, 116:2.04641E-04}
+    inAtomFractionsDict = {113:3.44262E-04, 115:7.68050E-03}
+    elemAtomFracDict = {'Ag':2.19403E-02, 'Cd':2.631098E-04, 'In':8.024762E-03}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, agAtomFractionsDict, 'Ag')
+    override_abundances(ZAList, abundanceDict, cdAtomFractionsDict, 'Cd')
+    override_abundances(ZAList, abundanceDict, inAtomFractionsDict, 'In')
+
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, atomDensity=atomDensity)
+    return material
+
+def get_CASL_fuel_p4_211_material():
+    shortName = 'CASLfuel_P4_211'
+    longName = 'CASL fuel P4 211 material'
+    massDensity = 10.257 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    uAtomFractionsDict = {234:4.04814E-06, 235:4.88801E-04, 236:2.23756E-06, 238:2.23844E-02}
+    elemAtomFracDict = {'O':4.57591E-02, 'U':2.28794867E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, uAtomFractionsDict, 'U')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_fuel_p4_262_material():
+    shortName = 'CASLfuel_P4_262'
+    longName = 'CASL fuel P4 262 material'
+    massDensity = 10.257 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    uAtomFractionsDict = {234:5.09503E-06, 235:6.06733E-04, 236:2.76809E-06, 238:2.22663E-02}
+    elemAtomFracDict = {'O':4.57617E-02, 'U':2.28808961E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, uAtomFractionsDict, 'U')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_moderator_p4_material():
+    shortName = 'CASLmod_P4'
+    longName = 'CASL moderator p4 material'
+    massDensity = 0.743 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'h2o'
+    bAtomFractionsDict = {10:1.12012E-05, 11:4.50862E-05}
+    elemAtomFracDict = {'H':4.96194E-02, 'O':2.48097E-02, 'B':5.62874E-05}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_topnozzle_material():
+    shortName = 'CASLtopnozzle'
+    longName = 'CASL top nozzle material'
+    atomDensity = 7.71078E-02 #atom/barn-cm
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    hAtomFractionsDict  = {1:4.01187E-02}
+    bAtomFractionsDict  = {10:9.05410E-06, 11:3.64439E-05}
+    siAtomFractionsDict = {28:3.02920E-04, 29:1.53886E-05, 30:1.01561E-05}
+    crAtomFractionsDict = {50:1.46468E-04, 52:2.82449E-03, 53:3.20275E-04, 54:7.97232E-05}
+    feAtomFractionsDict = {54:6.60188E-04, 56:1.03635E-02, 57:2.39339E-04, 58:3.18517E-05}
+    niAtomFractionsDict = {58:1.01650E-03, 60:3.91552E-04, 61:1.70205E-05, 62:5.42688E-05, 64:1.38207E-05}
+    elemAtomFracDict = {'H':4.01187E-02, 'B':4.5498E-05, 'C':6.14459E-05, 'O':2.00593E-02, 'Si':3.28465E-04, 'P':1.34026E-05, 'Cr':3.3709562E-03, 'Mn':3.35836E-04, 'Fe':1.12949E-02, 'Ni':1.47934E-03}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, hAtomFractionsDict, 'H')
+    override_abundances_as_elemental(ZAList, abundanceDict,  'C')
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    override_abundances(ZAList, abundanceDict, siAtomFractionsDict, 'Si')
+    override_abundances(ZAList, abundanceDict, crAtomFractionsDict, 'Cr')
+    override_abundances(ZAList, abundanceDict, feAtomFractionsDict, 'Fe')
+    override_abundances(ZAList, abundanceDict, niAtomFractionsDict, 'Ni')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=atomDensity)
+    return material
+
+def get_CASL_bottomnozzle_material():
+    shortName = 'CASLbottomnozzle'
+    longName = 'CASL bottom nozzle material'
+    atomDensity = 8.26822E-02 #atom/barn-cm
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    hAtomFractionsDict  = {1:3.57638E-02}
+    bAtomFractionsDict  = {10:8.07351E-06, 11:3.24969E-05}
+    siAtomFractionsDict = {28:4.41720E-04, 29:2.24397E-05, 30:1.48097E-05}
+    crAtomFractionsDict = {50:2.13581E-04, 52:4.11869E-03, 53:4.67027E-04, 54:1.16253E-04}
+    feAtomFractionsDict = {54:9.62690E-04, 56:1.51122E-02, 57:3.49006E-04, 58:4.64463E-05}
+    niAtomFractionsDict = {58:1.48226E-03, 60:5.70964E-04, 61:2.48194E-05, 62:7.91351E-05, 64:2.01534E-05}
+    elemAtomFracDict = {'H':4.01187E-02, 'B':4.05704E-05, 'C':8.96008E-05, 'O':1.78819E-02, 'Si':4.78969E-04, 'P':1.95438E-05, 'Cr':4.91555E-03, 'Mn':4.89719E-04, 'Fe':1.64703E-02, 'Ni':2.17733E-03}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, hAtomFractionsDict, 'H')
+    override_abundances_as_elemental(ZAList, abundanceDict,  'C')
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    override_abundances(ZAList, abundanceDict, siAtomFractionsDict, 'Si')
+    override_abundances(ZAList, abundanceDict, crAtomFractionsDict, 'Cr')
+    override_abundances(ZAList, abundanceDict, feAtomFractionsDict, 'Fe')
+    override_abundances(ZAList, abundanceDict, niAtomFractionsDict, 'Ni')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=atomDensity)
+    return material
+
+def get_CASL_coreplates_material():
+    shortName = 'CASLcoreplates'
+    longName = 'CASL core plates material'
+    atomDensity = 7.97787E-02 #atom/barn-cm
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    hAtomFractionsDict  = {1:2.48098E-02}
+    bAtomFractionsDict  = {10:5.62115E-06, 11:2.26258E-05}
+    siAtomFractionsDict = {28:7.90985E-04, 29:4.01826E-05, 30:2.65197E-05}
+    crAtomFractionsDict = {50:3.82458E-04, 52:7.37532E-03, 53:8.36302E-04, 54:2.08173E-04}
+    feAtomFractionsDict = {54:1.72388E-03, 56:2.70613E-02, 57:6.24963E-04, 58:8.31710E-05}
+    niAtomFractionsDict = {58:2.65427E-03, 60:1.02242E-03, 61:4.44439E-05, 62:1.41707E-04, 64:3.60885E-05}
+    elemAtomFracDict = {'H':2.48098E-02, 'B':2.82470E-05, 'C':1.60447E-04, 'O':1.24049E-02, 'Si':8.57687E-04, 'P':3.49969E-05, 'Cr':7.96595E-03, 'Mn':8.76936E-04, 'Fe':2.87852E-02, 'Ni':3.85449E-03}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, hAtomFractionsDict, 'H')
+    override_abundances_as_elemental(ZAList, abundanceDict,  'C')
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    override_abundances(ZAList, abundanceDict, siAtomFractionsDict, 'Si')
+    override_abundances(ZAList, abundanceDict, crAtomFractionsDict, 'Cr')
+    override_abundances(ZAList, abundanceDict, feAtomFractionsDict, 'Fe')
+    override_abundances(ZAList, abundanceDict, niAtomFractionsDict, 'Ni')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=atomDensity)
+    return material
+
+def get_CASL_fuel_p5_31_material():
+    shortName = 'CASLfuel_P5_31'
+    longName = 'CASL fuel P5 material'
+    massDensity = 10.257 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    uAtomFractionsDict = {234:6.11864E-06, 235:7.18132E-04, 236:03.29861E-06, 238:2.21546E-02}
+    elemAtomFracDict = {'O':4.57642E-02, 'U':2.28821E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, uAtomFractionsDict, 'U')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_moderator_p5_material():
+    shortName = 'CASLmod_P5'
+    longName = 'CASL moderator P5 material'
+    massDensity = 0.743 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'h2o'
+    bAtomFractionsDict = {10:1.06329E-05, 11:4.27988E-05}
+    elemAtomFracDict = {'H':4.96228E-02, 'O':2.48114E-02, 'B':5.34317E-05}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_pyrex_p5_material():
+    shortName = 'CASLpyrex_P5'
+    longName = 'CASL PYREX P5 material'
+    massDensity = 2.25 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    bAtomFractionsDict = {10:9.61468E-04, 11:3.89444E-03}
+    siAtomFractionsDict = {28:1.81641E-02, 29:9.22749E-04, 30:6.08994E-04}
+    elemAtomFracDict = {'B':4.85591E-03, 'O':4.66888E-02, 'Si':1.96958E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    override_abundances(ZAList, abundanceDict, siAtomFractionsDict, 'Si')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_StainlessSteel_p5_material():
+    shortName = 'CASLss_P5'
+    longName = 'CASL stainless-steel P5 material'
+    massDensity = 8.0 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    siAtomFractionsDict = {28:1.58197E-03, 29:8.03653E-05, 30:5.30394E-05}
+    crAtomFractionsDict = {50:7.64915E-04, 52:1.47506E-02, 53:1.67260E-03, 54:4.16346E-04}
+    feAtomFractionsDict = {54:3.44776E-03, 56:5.41225E-02, 57:1.24992E-03, 58:1.66342E-04}
+    niAtomFractionsDict = {58:5.30854E-03, 60:2.04484E-03, 61:8.88879E-05, 62:2.83413E-04, 64:7.21770E-05}
+    elemAtomFracDict = {'C':3.20895E-04, 'Si':1.71537E-03, 'P':6.99938E-05, 'Cr':1.76045E-02, 'Mn':1.75387E-03, 'Fe':5.89865E-02, 'Ni':7.79786E-03}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances_as_elemental(ZAList, abundanceDict,  'C')
+    override_abundances(ZAList, abundanceDict, siAtomFractionsDict, 'Si')
+    override_abundances(ZAList, abundanceDict, crAtomFractionsDict, 'Cr')
+    override_abundances(ZAList, abundanceDict, feAtomFractionsDict, 'Fe')
+    override_abundances(ZAList, abundanceDict, niAtomFractionsDict, 'Ni')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_CASL_B4C_p5_material():
+    shortName = 'CASLb4c_P5'
+    longName = 'CASL B4C P5 material'
+    atomDensity = 9.59100E-02 #1/barn-cm
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    bAtomFractionsDict = {10:1.52689E-02, 11:6.14591E-02}
+    elemAtomFracDict = {'B':7.67280E-02, 'C':1.91820E-02}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
+    override_abundances_as_elemental(ZAList, abundanceDict, 'C')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, atomDensity=atomDensity)
+    return material
+
+def get_CASL_AIC_p5_material():
+    shortName = 'CASLaic_P5'
+    longName = 'CASL AIC P5 material'
+    atomDensity = 5.38440718E-02 #1/barn-cm
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    agAtomFractionsDict = {107:2.36159E-02, 109:2.19403E-02}
+    cdAtomFractionsDict = {106:3.41523E-05, 108:2.43165E-05, 110:3.41250E-04, 111:3.49720E-04, 112:6.59276E-04, 113:3.33873E-04, 114:7.84957E-04, 116:2.04641E-04}
+    inAtomFractionsDict = {113:3.44262E-04, 115:7.68050E-03}
+    elemAtomFracDict = {'Ag':2.19403E-02, 'Cd':2.631098E-04, 'In':8.024762E-03}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, agAtomFractionsDict, 'Ag')
+    override_abundances(ZAList, abundanceDict, cdAtomFractionsDict, 'Cd')
+    override_abundances(ZAList, abundanceDict, inAtomFractionsDict, 'In')
+
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, atomDensity=atomDensity)
+    return material
+
+def get_CASL_gas_p5_material():
+    shortName = 'CASLgas_P5'
+    longName = 'CASL gas P5 material'
+    atomDensity = 2.68714e-05
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    elemAtomFracDict = {'He':1.0}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, atomDensity=atomDensity)
+    return material
+
+def get_CASL_cladding_p5_material():
+    shortName = 'CASLclad_P5'
+    longName = 'CASL cladding P5 material'
+    massDensity = 6.56 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 565. #K
+    thermalOpt = 'free'
+    crAtomFractionsDict = {50:3.30121E-06, 52:6.36606E-05, 53:7.21860E-06, 54:1.79686E-06}
+    feAtomFractionsDict = {54:8.68307E-06, 56:1.36306E-04, 57:3.14789E-06, 58:4.18926E-07}
+    zrAtomFractionsDict = {90:2.18865E-02, 91:4.77292E-03, 92:7.29551E-03, 94:7.39335E-03, 
+                           96:1.19110E-03}
+    snAtomFractionsDict = {112:4.68066E-06, 114:3.18478E-06, 115:1.64064E-06, 116:7.01616E-05,
+                           117:3.70592E-05, 118:1.16872E-04, 119:4.14504E-05, 120:1.57212E-04,
+                           122:2.23417E-05, 124:2.79392E-05}
+    hfAtomFractionsDict = {174:3.54138E-09, 176:1.16423E-07, 177:4.11686E-07, 178:6.03806E-07,
+                           179:3.01460E-07, 180:7.76449E-07}
+    elemAtomFracDict = {'Cr':7.59773E-05, 'Fe':1.48556E-04, 'Zr':4.25394E-02, 'Sn':4.82542E-04, 'Hf':2.21337E-06}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, crAtomFractionsDict, 'Cr')
+    override_abundances(ZAList, abundanceDict, feAtomFractionsDict, 'Fe')
+    override_abundances(ZAList, abundanceDict, zrAtomFractionsDict, 'Zr')
+    override_abundances(ZAList, abundanceDict, snAtomFractionsDict, 'Sn')
+    override_abundances(ZAList, abundanceDict, hfAtomFractionsDict, 'Hf')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
+def get_410Design_material():
+    shortName = '410Design'
+    longName = '410 Design Project'
+    massDensity = 10.257 #g/cc
+    fuelRadius = 0.4096 #cm
+    temperature = 600. #K
+    thermalOpt = 'h2o'
+    uAtomFractionsDict = {235:7.18132E-04, 238:03.29861E-06, 239:2.21546E-02}
+    npAtomFractionsDict = {239:1.0}
+    puAtomFractionsDict = {239:1.0}
+    elemAtomFracDict = {'O':4.57642E-02, 'H': 9.15E-02, 'U':2.28821E-02, 'Np': 2.2E-04, 'Pu': 2.2E-03}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, uAtomFractionsDict, 'U')
+    override_abundances(ZAList, abundanceDict, npAtomFractionsDict, 'Np')
+    override_abundances(ZAList, abundanceDict, puAtomFractionsDict, 'Pu')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, massDensity=massDensity)
+    return material
+
 def get_cert_material():
     shortName = 'cert'
     longName = 'CERT material'
@@ -344,7 +927,7 @@ def get_cert_Hpoly_Cfree():
     return material
 
 def get_cert_Hh2o():
-    shortName = 'certHhwo'
+    shortName = 'certHh2o'
     longName = 'CERT H-h2o'
     massDensity = 5. #g/cc
     fuelRadius = 0. #cm
@@ -1096,6 +1679,40 @@ def get_triga_fuel_material():
         temperatureIndex=temperatureIndex)
     return material
 
+def get_ctriga_fuel_material():
+    shortName = 'tcFUEL'
+    longName = 'U-ZrH fuel complete'
+    atomDensity = 8.71115E-2
+    fuelRadius = 1.7411 - 0.3175 #cm
+    temperature = 296. #K
+    temperatureIndex = 0 # X in .9Xc
+    thermalOpt = 'zrh'
+    uAtomFractionsDict = {234: 9.44763913e-05, 235: 1.24205225e-02 , 236: 1.38902530e-04, 237: 7.14717721e-25, 238: 4.96062106e-02}
+    hfAtomFractionsDict = {176: 1.16946632e-06, 177: 4.13537520e-06, 178: 6.06521696e-06, 179: 3.02816184e-06, 180: 7.81190561e-06}
+    erAtomFractionsDict = {164: 4.24355777e-05, 166: 8.89255527e-04, 167:  6.06158167e-04, 168: 7.15069965e-04, 170: 3.95199540e-04}
+    zrAtomFractionsDict = {90: 1.90650598e-01, 91: 4.15762821e-02, 92: 6.35501995e-02, 94: 6.44024762e-02, 96: 1.03755428e-02}
+    hAtomFractionsDict = {1: 5.64241953e-01, 2: 6.48952875e-05}
+    elemAtomFracDict = {'U': 6.22603E-02, 'Zr': 3.70556E-01, 'Hf': 2.22102E-05, 'Er': 2.64813E-03, 'C': 2.05142E-04, 'H': 5.64308E-01}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances(ZAList, abundanceDict, uAtomFractionsDict, 'U')
+    override_abundances(ZAList, abundanceDict, hfAtomFractionsDict, 'Hf')
+    override_abundances(ZAList, abundanceDict, erAtomFractionsDict, 'Er')
+    override_abundances(ZAList, abundanceDict, zrAtomFractionsDict, 'Zr')
+    override_abundances(ZAList, abundanceDict, hAtomFractionsDict, 'H')
+    override_abundances_as_elemental(ZAList, abundanceDict, 'C')
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemAtomFracDict=elemAtomFracDict, atomDensity=atomDensity,
+        temperatureIndex=temperatureIndex)
+    return material
+
 def get_triga_clad_material():
     shortName = 'tCLAD'
     longName = 'SS304 clad'
@@ -1119,10 +1736,45 @@ def get_triga_clad_material():
         temperatureIndex=temperatureIndex)
     return material
 
+def get_ctriga_clad_material():
+    shortName = 'tcCLAD'
+    longName = 'SS304 clad complete'
+    massDensity = 8.0
+    fuelRadius = 1.7920 #cm
+    temperature = 296. #K
+    temperatureIndex = 0 # X in .9Xc
+    thermalOpt = 'free'
+    siAtomFractionsDict = {28: 0.004594, 29: 0.000241, 30: 0.000165}
+    sAtomFractionsDict  = {32: 0.000142, 33: 0.000001, 34: 0.000007}
+    crAtomFractionsDict = {50: 0.007930, 52: 0.159031, 53: 0.018378, 54: 0.004661}
+    feAtomFractionsDict = {54: 0.039996, 56: 0.644764, 57: 0.015026, 58: 0.002039}
+    niAtomFractionsDict = {58: 0.062340, 60: 0.024654, 60: 0.001085, 62: 0.003504, 64: 0.000917}
+    elemMassFracDict = {'C': 0.0003, 'Si': 0.005, 'P': 0.000225, 'S': 0.00015, 'Cr': 0.19, 'Mn': 0.01, 'Fe': 0.701825, 'Ni': 0.0925}
+    #
+    chordLength = calc_chord_length(fuelRadius)
+    symDict, ZList, ZAList = get_all_isotopes(elemMassFracDict)
+    abundanceDict = lookup_natl_abundances(ZAList)
+    override_abundances_as_elemental(ZAList, abundanceDict, 'C')
+    override_abundances(ZAList, abundanceDict, siAtomFractionsDict, 'Si', 'Mass')
+    override_abundances(ZAList, abundanceDict, sAtomFractionsDict, 'S', 'Mass')
+    override_abundances(ZAList, abundanceDict, crAtomFractionsDict, 'Cr', 'Mass')
+    override_abundances(ZAList, abundanceDict, feAtomFractionsDict, 'Fe', 'Mass')
+    override_abundances(ZAList, abundanceDict, niAtomFractionsDict, 'Ni', 'Mass')
+
+    #
+    material = Material(
+        shortName=shortName, longName=longName,
+        temperature=temperature, thermalOpt=thermalOpt,
+        symDict=symDict, ZList=ZList, ZAList=ZAList,
+        abundanceDict=abundanceDict, chordLength=chordLength,
+        elemMassFracDict=elemMassFracDict, massDensity=massDensity,
+        temperatureIndex=temperatureIndex)
+    return material
+
 def get_triga_zirconium_material():
     shortName = 'tZIRC'
     longName = 'Zr fuel center'
-    atomDensity = 4.29600E-2
+    atomDensity =  0.01296
     fuelRadius = 0.3175 #cm
     temperature = 296. #K
     temperatureIndex = 0 # X in .9Xc
@@ -1198,12 +1850,14 @@ def get_triga_b4c_material():
     temperature = 296. #K
     temperatureIndex = 0 # X in .9Xc
     thermalOpt = 'free'
-    elemAtomFracDict = {'C': 0.201616066, 'B': 0.798383934}
+    bAtomFractionsDict = {10: 0.020950, 11: 0.084310}
+    elemAtomFracDict = {'C': 0.026320, 'B': 0.105260}
     #
     chordLength = calc_chord_length(fuelRadius)
     symDict, ZList, ZAList = get_all_isotopes(elemAtomFracDict)
     abundanceDict = lookup_natl_abundances(ZAList)
     override_abundances_as_elemental(ZAList, abundanceDict, 'C')
+    override_abundances(ZAList, abundanceDict, bAtomFractionsDict, 'B')
     #
     material = Material(
         shortName=shortName, longName=longName,
@@ -1240,7 +1894,7 @@ def get_triga_irradiation_tube_material():
 def get_triga_moderator_material():
     shortName = 'tMOD'
     longName = 'light water'
-    atomDensity = 1.00037E-1
+    atomDensity = 0.09989 #1.00037E-1
     fuelRadius = 1.7920 #cm
     temperature = 293.6 #K
     temperatureIndex = 0 # X in .9Xc
@@ -1576,6 +2230,7 @@ def get_multi_temperature_h2o_material_base():
     fuelRadius = 1.7920 #cm
     temperature = 293.6 #K
     thermalOpt = 'h2o'
+    temperatureIndex = 0 # X in .9Xc
     # Hack to align thermal grids:
     #thermalOpt = 'free'
     elemAtomFracDict = {'H': 2, 'O': 1}
@@ -1777,6 +2432,37 @@ def get_multi_temperature_depleted_triga_fuel_material_T7():
 
 ###############################################################################
 
+def get_multi_temperature_triga_b4c_material_base():
+    return get_triga_b4c_material()
+
+def get_multi_temperature_triga_b4c_material_Tgrid():
+    return [296, 400, 500, 600, 700, 800, 1000, 1200] #K
+
+def get_multi_temperature_triga_b4c_material_Ti(iT):
+    material = get_multi_temperature_triga_b4c_material_base()
+    Tgrid = get_multi_temperature_triga_b4c_material_Tgrid()
+    T = Tgrid[iT]
+    shortName = 'mtTB4C_{}'.format(iT)
+    longName = '{} ({} K)'.format(material.longName, T)
+    material.update_temperature(T)
+    material.update_temperature_index(iT) # X in .9Xc
+    # MASS DENSITY SHOULD BE UPDATED TO ACCOUNT FOR THERMAL EXPANSION
+    #material.update_mass_density(rho)
+    material.update_names(shortName, longName)
+    return material
+
+def get_multi_temperature_triga_b4c_material_T0():
+    return get_multi_temperature_triga_b4c_material_Ti(0)
+
+def get_multi_temperature_triga_b4c_material_T5():
+    return get_multi_temperature_triga_b4c_material_Ti(5)
+
+def get_multi_temperature_triga_b4c_material_T7():
+    return get_multi_temperature_triga_b4c_material_Ti(7)
+
+
+###############################################################################
+
 def get_all_isotopes(elemDict):
     symList = elemDict.keys()
     symDict = {}
diff --git a/src/write_njoy.py b/src/write_njoy.py
index d4d5e30..5d08eb7 100644
--- a/src/write_njoy.py
+++ b/src/write_njoy.py
@@ -23,7 +23,7 @@ def create_thermal_ace_njoy_script(dat, tapes):
     endfDirr = dirDict['endf']
     pendfRootDirr = dirDict['pendf']
     aceRootDirr = dirDict['ace']
-    njoyPath = os.path.join(dirDict['njoyInstall'], 'xnjoy')
+    njoyPath = os.path.join(dirDict['njoyInstall'], 'xnjoy_thermr') #using the latest patched NJOY2012
     #
     pendfDirr = os.path.join(pendfRootDirr, str(dat.nuclideName))
     aceDirr = os.path.join(aceRootDirr, str(dat.thermalFileName))
@@ -79,6 +79,7 @@ def create_thermal_ace_njoy_script(dat, tapes):
         deck.append(["./xnjoy<input"])
     deck.append(["echo 'Cleaning up and saving ACE files'"])
     aceFiles = []
+    aceCERTDir = '/scratch/yunhuang/barnfire/xs/cert_ace'
     for i in range(numTemperatures):
         tapeAceOut, tapeXSDirOut, viewrOut = tapes.aceStart + i, tapes.aceXSDirStart + i, tapes.viewrAceStart + i
         aceFile = aceFileTemplate.format(H=dat.thermalFileName, s=i, e=dat.aceExt)
@@ -87,6 +88,8 @@ def create_thermal_ace_njoy_script(dat, tapes):
         plotName = 'p_xs_{}.pdf'.format(aceFile)
         deck.append(['cp -f tape{0:02g} {1}'.format(tapeAceOut, aceFile)])
         deck.append(['cp -f tape{0:02g} {1}'.format(tapeXSDirOut, aceFileDir)])
+        deck.append(['cp -f {0} {1}/'.format(aceFile, aceCERTDir)])
+        deck.append(['cp -f {0} {1}/'.format(aceFileDir, aceCERTDir)])
         deck.append(['ps2pdf tape{0:02g} {1}'.format(viewrOut, plotName)])
     deck.append(['rm -f xnjoy'])
     tapeNamesToRemove = ' '.join(get_temporary_tapes_ace(tapes, pendfOutTape, endfNonFreeThermalTapeList, numTemperatures, False))
@@ -106,7 +109,7 @@ def create_njoy_script(dat, tapes):
     pendfRootDirr = dirDict['pendf']
     gendfRootDirr = dirDict['gendf']
     aceRootDirr = dirDict['ace']
-    njoyPath = os.path.join(dirDict['njoyInstall'], 'xnjoy')
+    njoyPath = os.path.join(dirDict['njoyInstall'], 'xnjoy_thermr')
     #
     pendfDirr = os.path.join(pendfRootDirr, str(dat.nuclideName))
     gendfDirr = os.path.join(gendfRootDirr, str(dat.nuclideName))
@@ -145,6 +148,9 @@ def create_njoy_script(dat, tapes):
     relPathPendfGendf = os.path.relpath(pendfPath, gendfDirr)
     readgrouprPath = os.path.join(pyDirr, 'Readgroupr.py')
     #
+    matxsFile = 'matxs' 
+    matxsDir = '/scratch/yunhuang/barnfire/xs/cert_matxs'
+    #
     aceScriptFile = 'runNJOY.sh'
     aceScriptPath = os.path.join(aceDirr, aceScriptFile)
     aceFileTemplate = '{Z}{A:03}.{e}{s}c'
@@ -185,7 +191,7 @@ def create_njoy_script(dat, tapes):
     try_mkdir(pendfDirr)
     print_njoy_file(pendfScriptPath, deck)
 
-    # Script to run NJOY from PENDF to GENDF:
+    # Script to run NJOY from PENDF to MATXS:
     doPlotr = False
     deck = []
     deck.append(["#! /usr/bin/env bash"])
@@ -199,26 +205,31 @@ def create_njoy_script(dat, tapes):
     create_moder_input(deck, dat, tapes.endf, tapes.bendf)
     dat.hasDelayedNu, dat.hasDelayedChis = hasDelayed(endfPath)
     create_groupr_input(deck, dat, tapes.bendf, pendfOutTape, tapes.grouprIn, tapes.grouprOut)
-    create_moder_input(deck, dat, tapes.grouprOut, tapes.gendfOut)
+    create_moder_input(deck, dat, tapes.grouprOut, tapes.gendfOut)  # get GENDF in ASCII for record
+    create_matxsr_input(deck, dat, tapes.gendfOut, tapes.gaminrOut, tapes.matxsOut)
+    # It seems MODER is not needed after MATXSR
+    #  create_moder_input(deck, dat, tapes.matxsOut, tapes.gendfOut)
     if False and dat.numGroups <= 200:
         doPlotr = True
         plotNames = create_plotr_inputs(deck, dat, tapes.gendfOut, tapes.plotrOut, tapes.viewrOut)
     deck.append(['stop'])
     deck.append(["EOF"])
     deck.append(["./xnjoy<input"])
-    deck.append(["echo 'Cleaning up and saving GROUPR file'"])
+    deck.append(["echo 'Cleaning up, saving GROUPR and MATXS file'"])
     if doPlotr:
         for plotName, viewrOut in zip(plotNames, tapes.viewrOut):
             deck.append(['ps2pdf tape{0:02g} {1}'.format(viewrOut, plotName)])
     deck.append(['cp -f tape{0:02g} {1}'.format(tapes.gendfOut, gendfFile)])
+    deck.append(['cp -f tape{0:02g} {1}'.format(tapes.matxsOut, matxsFile)])
+    deck.append(['cp -f matxs {0}/{1}.dat'.format(matxsDir, dat.nuclideName)])
     deck.append(['rm -f xnjoy'])
     tapeNamesToRemove = ' '.join(get_temporary_tapes_gendf(tapes, pendfOutTape, doPlotr))
     deck.append(['rm -f {0}'.format(tapeNamesToRemove)])
     optionalComment = ''
     if True and dat.numGroups > 200:
         optionalComment = '#'
-    deck.append(["{0}echo 'Reading in XS file and pickling for future use'".format(optionalComment)])
-    deck.append(['{0} {1} -I {2} -p all -f csc'.format(optionalComment, readgrouprPath, gendfFile)])
+ #   deck.append(["{0}echo 'Reading in XS file and pickling for future use'".format(optionalComment)])
+ #   deck.append(['{0} {1} -I {2} -p all -f csc'.format(optionalComment, readgrouprPath, gendfFile)])
     try_mkdir(gendfDirr)
     print_njoy_file(gendfScriptPath, deck)
 
@@ -252,6 +263,7 @@ def create_njoy_script(dat, tapes):
         deck.append(["./xnjoy<input"])
     deck.append(["echo 'Cleaning up and saving ACE files'"])
     aceFiles = []
+    aceCERTDir = '/scratch/yunhuang/barnfire/xs/cert_ace'
     for i in range(numTemperatures):
         tapeAceOut, tapeXSDirOut, viewrOut = tapes.aceStart + i, tapes.aceXSDirStart + i, tapes.viewrAceStart + i
         aceFile = aceFileTemplate.format(A=dat.A, Z=dat.Z, s=i, e=dat.aceExt)
@@ -260,6 +272,8 @@ def create_njoy_script(dat, tapes):
         plotName = 'p_xs_{}.pdf'.format(aceFile)
         deck.append(['cp -f tape{0:02g} {1}'.format(tapeAceOut, aceFile)])
         deck.append(['cp -f tape{0:02g} {1}'.format(tapeXSDirOut, aceFileDir)])
+        deck.append(['cp -f {0} {1}/'.format(aceFile, aceCERTDir)])
+        deck.append(['cp -f {0} {1}/'.format(aceFileDir, aceCERTDir)])
         deck.append(['ps2pdf tape{0:02g} {1}'.format(viewrOut, plotName)])
     deck.append(['rm -f xnjoy'])
     tapeNamesToRemove = ' '.join(get_temporary_tapes_ace(tapes, pendfOutTape, [], numTemperatures, True))
@@ -357,7 +371,7 @@ def create_reconr_input(deck, dat, tapeENDFIn, tapePENDFOut):
     deck.append((tapeENDFIn, tapePENDFOut))
     deck.append(("'pendf tape for " + dat.nuclideName + " from " + dat.endfName + "'",'/'))
     deck.append((dat.mat, ncards, '/'))
-    deck.append((dat.errTol, 0, dat.errTolMax, '/'))
+    deck.append((dat.errTol, '/'))
     deck.append(("'" + dat.nuclideName + " from " + dat.endfName + " tape'",'/'))
     deck.append(("'processed by the njoy nuclear data processing system'",'/'))
     deck.append(("'see original " + dat.endfName + " tape for details of evaluation'",'/'))
@@ -409,7 +423,7 @@ def create_thermr_input(deck, dat, index, tapeENDFThermalIn, tapePENDFIn, tapePE
     # For NJOY 2012, before numAtom, insert iform (probably should be 0).
     # NB NJOY2012 messes up the output of the thermal data
     # in the GENDF file, which messes up Readgroupr.py
-    deck.append((matThermal, dat.mat, dat.numAngles, len(dat.thermList), inelasticOpt, elasticOpt, numAtom, inelasticMT, 1))
+    deck.append((matThermal, dat.mat, dat.numAngles, len(dat.thermList), inelasticOpt, elasticOpt, dat.inelasForm, numAtom, inelasticMT, 1))
     deck.append(dat.thermList)
     deck.append((dat.thermrTol, dat.Emax))
 
@@ -546,6 +560,7 @@ def create_groupr_input(deck, dat, tapeENDFIn, tapePENDFIn, tapeGroupsIn, tapeGE
     namedRxts = []
     #
     namedRxts.append((3, '/'))
+    # inverse velocity not used for CERT
     namedRxts.append((3, 259, "'invel'", '/'))
     for mt in dat.thermalMTList:
         namedRxts.append((3,mt,"'(thermal)'",'/'))
@@ -564,7 +579,8 @@ def create_groupr_input(deck, dat, tapeENDFIn, tapePENDFIn, tapeGroupsIn, tapeGE
         # To get chi, currently need (6,18), so use this even if not including other MF6
         namedRxts.append((6, 18,"'fission'",'/'))
 
-    namedRxts = sorted(namedRxts)
+    # Sorted namedRxts will cause random string error in MATXS
+    ## namedRxts = sorted(namedRxts)
     # repeat total, elastic, fission, capture, thermal for each temperature
     toRepeat = [(3,1,"'total'",'/'), (3,2,"'elastic scat'",'/'), (3,102,"'rad capture'",'/')]
     if dat.isFissionable:
@@ -609,6 +625,8 @@ def create_groupr_input(deck, dat, tapeENDFIn, tapePENDFIn, tapeGroupsIn, tapeGE
         groupFormat.append('/')
         if rowPos != 1:
             deck.append(groupFormat)
+    if (dat.weightOpt == 4):
+        deck.append([0.4, 0.0255, 2.0e7, 8.2030E+05])
     for rxt in namedRxts:
         deck.append(rxt)
     deck.append((0, '/'))
@@ -618,6 +636,25 @@ def create_groupr_input(deck, dat, tapeENDFIn, tapePENDFIn, tapeGroupsIn, tapeGE
         deck.append((0, '/'))
     deck.append((0, '/'))
 
+def create_matxsr_input(deck, dat, tapeGENDFIn, tapeGAMINRIn, tapeMATXSOut):
+    groupMap = [-1, -1, 240, 30, 27, 50, 68, 100, 35, 69, 187, 70, 620, 80, 100, 640, 174, 175]
+    dat.numGroups = groupMap[dat.groupOpt]
+    if dat.groupOpt == 1:
+        dat.numGroups = len(dat.groupBdrs) - 1
+
+    deck.append(['matxsr'])
+    deck.append((tapeGENDFIn, tapeGAMINRIn, tapeMATXSOut, '/'))
+    deck.append((1, "'y.zhang cert data'", '/'))
+    deck.append((1, 2, 1, 1))
+    deck.append(["'" + str(dat.numGroups) + '-group ' + dat.nuclideName + ' neutron lib' + "'"])
+    deck.append(["'n'"])
+    deck.append([str(dat.numGroups)])
+    deck.append(("'nscat'", "'ntherm'", '/'))
+    deck.append((1, 1))
+    deck.append((1, 1))
+    deck.append(("'" + dat.nuclideName + "'", dat.mat, dat.mat))
+    
+
 ###############################################################################
 def print_njoy_file(filePath, deck):
     """Given a file path and deck, print out an NJOY file."""
@@ -692,6 +729,8 @@ def __init__(self):
         self.plotrOut = [31, 33, 35, 37]
         self.viewrOut = [32, 34, 36, 38]
         self.purrOut = -39
+        self.matxsOut = 42
+        self.gaminrOut = 0  #no photon data for CERT project
 
 
 def get_temporary_tapes_pendf(tapes, endfThermalTapeList, numThermals):
@@ -753,7 +792,7 @@ def get_temporary_tapes_ace(tapes, pendfOutTape, endfThermalTapeList, numTempera
 
 ###############################################################################
 class NJOYDat():
-    def __init__(self, A=None, Z=None, nuclideName=None, thermalNuclideName=None, thermalFileName=None, endfName=None, endfFile=None, mat=None, aceExt=None, thermList=[0], sig0List=[1e10], groupBdrs=None, groupOpt=0, isFissionable=False, hasDelayedNu=False, hasDelayedChis=False, includeMF6=True, usePURR=True, legendreOrder=0, inelasticMTList=[], thermalMTList=[], matThermalList=[], endfThermalFileList=[]):
+    def __init__(self, A=None, Z=None, nuclideName=None, thermalNuclideName=None, thermalFileName=None, endfName=None, endfFile=None, mat=None, aceExt=None, thermList=[0], sig0List=[1e10], groupBdrs=None, groupOpt=0, isFissionable=False, hasDelayedNu=False, hasDelayedChis=False, includeMF6=True, usePURR=False, legendreOrder=0, inelasticMTList=[], thermalMTList=[], matThermalList=[], endfThermalFileList=[]):
         # General
         # Nuclide information
         self.A = A
@@ -786,11 +825,13 @@ def __init__(self, A=None, Z=None, nuclideName=None, thermalNuclideName=None, th
         self.thermalMTList = thermalMTList
         # tolerance
         self.thermrTol = 0.001
+        # inelastic distribution form
+        self.inelasForm = 0
         # maximum energy for thermal treatment (eV) (originally used 1.9; others use 4.0)
         # self.Emax = 2.9
         #self.Emax = 3.3 # works with wiggle room for SHEM-361 and edits-12
         #self.Emax = 3.9 # does not work for ""
-        self.Emax = 3.6 # works for ""
+        self.Emax = 100 # xnjoy_thermr patch allows emax up to 100eV
         #
         # GROUPR
         self.includeMF6 = includeMF6
@@ -801,6 +842,7 @@ def __init__(self, A=None, Z=None, nuclideName=None, thermalNuclideName=None, th
         # ign: 1 for custom structure (see manual)
         self.groupOpt = groupOpt
         # iwt: weighting spectrum (see manual)
+        # use iwt = 4 for CERT
         self.weightOpt = 5
         # How many Legendre orders of transfer matrices to output; n in Pn
         self.legendreOrder = legendreOrder
@@ -816,7 +858,11 @@ def __init__(self, A=None, Z=None, nuclideName=None, thermalNuclideName=None, th
         # Thermal ACE
         self.numMix = 1 #True except for BeO or C6H6 (Benzene)
         self.numThermalAceBins = 128
-        self.emaxThermalAce = 1000.
+        # self.emaxThermalAce = 1000.
+        # NJOY default to 1000eV, which provides more fidelity in thermal treatment
+        # However, with the latest xnjoy_thermr patch, emax > 20eV will result in
+        # segfault in MCNP runtime.
+        self.emaxThermalAce = 20.
         # ZAIDs for which the thermal ACE treatment applies
         self.ZAIDs = []
 
@@ -841,7 +887,7 @@ def get_num_atoms_in_molecule(inelasticMT):
     return numAtomsInMolecule[inelasticMT]
 
 def get_elastic_option(inelasticMT):
-    if inelasticMT == 221:
+    if inelasticMT == 221 or inelasticMT == 222:  # water (MT222) has not elastic component in NJOY2012
         # Free gas
         return 0
     else:
@@ -854,7 +900,7 @@ def get_inelastic_option(inelasticMT):
         return 1
     else:
         # Read S(alpha, beta)
-        return 4
+        return 2
     #else:
     #    return 0