diff --git a/mg-criticality/make_library.py b/mg-criticality/make_library.py
index 2c61700..ad3d25b 100755
--- a/mg-criticality/make_library.py
+++ b/mg-criticality/make_library.py
@@ -16,6 +16,27 @@ GROUP_STRUCT = {1: openmc.mgxs.EnergyGroups(group_edges=[0.0, 20.0e6]),
 ###############################################################################
 
 
+def check_it(scatter_matrix):
+    import scipy.special as ss
+    groups = scatter_matrix.shape[0]
+    orders = scatter_matrix.shape[2]
+    Nmu = 50
+    mu = np.linspace(-1, 1, Nmu)
+    f = np.zeros(shape=(groups, groups, Nmu))
+    neg = False
+    for gin in range(groups):
+        for gout in range(groups):
+            data = np.array(scatter_matrix[gin, gout, :], copy=True)
+            if np.sum(data) > 0:
+                data /= data[0]
+                for l in range(orders):
+                    f[gin, gout, :] += ((float(l) + 0.5) *
+                                        ss.eval_legendre(l, mu) * data[l])
+                if np.any(f[gin, gout, :] < 0):
+                    neg = True
+    return neg
+
+
 def plot_it(scatter_matrix):
     import matplotlib.pyplot as plt
     import scipy.special as ss
@@ -44,6 +65,9 @@ def plot_it(scatter_matrix):
 
 
 def set_it(name, groups, order, fission, nu, absorption, scatt, total, chi):
+    neg = check_it(scatt)
+    if neg:
+        print(name + ' is negative!')
     xsd = openmc.XSdata(name, groups)
     xsd.order = order
     if fission is not None:
@@ -305,7 +329,6 @@ def create_2g():
     scatter = np.rollaxis(scatter, 0, 3)
     total = [0.2208, 0.3360]
     chi = [0.575, 0.425]
-
     Pu = set_it('Pu', groups, 0, fiss, nu, absorption, scatter, total, chi)
     mg_cross_sections_file.add_xsdata(Pu)
 
@@ -462,7 +485,7 @@ def create_2g():
     capture = [0.0087078, 0.02518]
     absorption = np.add(capture, fiss)
     scatter = np.array(
-        [[[0.31980, 0.0045552],
+        [[[0.31980, 0.004555],
           [0.000000, 0.42410]]])
     scatter = np.rollaxis(scatter, 0, 3)
     total = [0.33588, 0.54628]
@@ -491,8 +514,8 @@ def create_2g():
     ###########################################################################
     # 5.2.2 UD2O-1
     nu = [2.50, 2.50]
-    fiss = np.array([0.0028172, 0.097])
-    capture = [0.0087078, 0.02518]
+    fiss = np.array([0.002817, 0.097])
+    capture = [0.008708, 0.02518]
     absorption = np.add(capture, fiss)
     scatter = np.array(
         [[[0.31980, 0.06694], [0.004555, -0.0003972]],
@@ -526,7 +549,7 @@ def create_3g():
           [0.000, 0.000, 2.000]]])
     scatter = np.rollaxis(scatter, 0, 3)
     total = [0.240, 0.975, 3.10]
-    chi = [0.96, 0.05, 0.0]
+    chi = [0.96, 0.04, 0.0]
 
     URR = set_it('URR', groups, 0, fiss, nu, absorption, scatter, total, chi)
     mg_cross_sections_file.add_xsdata(URR)
@@ -552,10 +575,10 @@ def create_6g():
           [0.000, 0.000, 2.000, 0.00, 0.00, 0.00],
           [0.000, 0.000, 0.000, 2.00, 0.00, 0.00],
           [0.000, 0.000, 0.000, 0.275, 0.60, 0.00],
-          [0.000, 0.000, 2.000, 0.033, 0.171, 0.024]]])
+          [0.000, 0.000, 0.000, 0.033, 0.171, 0.024]]])
     scatter = np.rollaxis(scatter, 0, 3)
     total = [0.240, 0.975, 3.10, 3.10, 0.975, 0.240]
-    chi = [0.48, 0.02, 0.0, 0.0, 0.02, 0.048]
+    chi = [0.48, 0.02, 0.0, 0.0, 0.02, 0.48]
 
     URR = set_it('URR', groups, 0, fiss, nu, absorption, scatter, total, chi)
     mg_cross_sections_file.add_xsdata(URR)
diff --git a/mg-criticality/make_model.py b/mg-criticality/make_model.py
index a80ebd7..44a249f 100644
--- a/mg-criticality/make_model.py
+++ b/mg-criticality/make_model.py
@@ -1,6 +1,5 @@
-#!/usr/bin/env python
-
 import numpy as np
+from subprocess import CalledProcessError
 
 import openmc
 import openmc.checkvalue as cv
@@ -55,7 +54,7 @@ class Case(object):
             mats.append(openmc.Material(name=self.mat_names[i]))
             mats[-1].set_density('macro', 1.0)
             mats[-1].add_macroscopic(macros[-1])
-            materials_file.add_material(mats[-1])
+            materials_file.append(mats[-1])
 
         materials_file.cross_sections = GROUP_FILES[self.groups]
 
@@ -101,6 +100,8 @@ class Case(object):
 
         settings_file.source = openmc.source.Source(space=uniform_dist)
 
+        settings_file.output = {'summary': False}
+
         return settings_file
 
     def make_tallies(self, r=None):
@@ -223,10 +224,10 @@ class Case(object):
             surfs = []
             for r, rad in enumerate(self.rad):
                 if r == len(self.rad) - 1:
-                    surfs.append(openmc.ZCylinder(R=rad,
+                    surfs.append(openmc.ZCylinder(r=rad,
                                                   boundary_type='vacuum'))
                 else:
-                    surfs.append(openmc.ZCylinder(R=rad))
+                    surfs.append(openmc.ZCylinder(r=rad))
 
             # Instantiate Cells
             cells = []
@@ -242,9 +243,9 @@ class Case(object):
             surfs = []
             for r, rad in enumerate(self.rad):
                 if r == len(self.rad) - 1:
-                    surfs.append(openmc.Sphere(R=rad, boundary_type='vacuum'))
+                    surfs.append(openmc.Sphere(r=rad, boundary_type='vacuum'))
                 else:
-                    surfs.append(openmc.Sphere(R=rad))
+                    surfs.append(openmc.Sphere(r=rad))
 
             # Instantiate Cells
             cells = []
@@ -301,9 +302,15 @@ class Case(object):
         geometry.export_to_xml()
 
     def execute(self, quiet=True):
-        returncode = openmc.run(output=(not quiet))
-        spfile = 'statepoint.' + str(self.batches) + '.h5'
-        sp = openmc.StatePoint(spfile)
-        self.keff = sp.k_combined[:]
+        success = True
+        try:
+            openmc.run(output=(not quiet))
+        except CalledProcessError:
+            success = False
 
-        return returncode
+        if success:
+            spfile = 'statepoint.' + str(self.batches) + '.h5'
+            sp = openmc.StatePoint(spfile, autolink=False)
+            self.keff = sp.k_combined
+
+        return success
diff --git a/mg-criticality/repo.py b/mg-criticality/repo.py
index 7dc2230..4dca3ec 100644
--- a/mg-criticality/repo.py
+++ b/mg-criticality/repo.py
@@ -1,12 +1,10 @@
-#!/usr/bin/env python
-
 from make_model import Case
 
-# GLOBAL DATA
-batches = 5000
-inactive = 500
-particles = 10000
-tab_leg = {'enable': True, 'num_points': 33}
+# GLOBAL OPTIONS
+batches = 1000
+inactive = 50
+particles = 2000
+tab_leg = {'enable': True}
 
 
 def build_cases():
@@ -574,22 +572,6 @@ def build_cases():
                       params))
     names[name] = len(cases) - 1
 
-    # # CASE -11
-    # case = -11
-    # name = 'Ua-1-1-IN'
-    # mat_names = ['Ua1']
-    # groups = 1
-    # order = 1
-    # geom = 'IN'
-    # rad = [1.]
-    # ref_k = 2.25
-    # mesh_dim = [20, 1, 1]
-    # params = {'mesh_dim': mesh_dim, 'tab_leg': tab_leg, 'batches': batches,
-    #           'inactive': inactive, 'particles': particles}
-    # cases.append(Case(case, name, mat_names, groups, order, geom, rad, ref_k,
-    #                   params))
-    # names[name] = len(cases) - 1
-
     # CASE 36
     case = 36
     name = 'Ua-1-1-CY'
diff --git a/mg-criticality/run_cases.py b/mg-criticality/run_cases.py
index d228d46..c346c89 100755
--- a/mg-criticality/run_cases.py
+++ b/mg-criticality/run_cases.py
@@ -1,20 +1,18 @@
 #!/usr/bin/env python
 
-from __future__ import print_function
 from optparse import OptionParser
-
-import matplotlib.pyplot as plt
-import numpy as np
+import csv
 
 import repo
 import make_library
 
 
 def print_case(case):
-    bias = np.abs(1.0E5 * (case.keff[0] - case.ref_k))
-    print('\tCalculated keff = {0:1.6f}'.format(case.keff[0]))
+    bias = 1.0E5 * (case.keff - case.ref_k)
+    print('\tCalculated keff = {0:1.6f}'.format(case.keff))
     print('\tReference keff  = {0:1.6f}'.format(case.ref_k))
     print('\tBias [pcm]      = {0:1.1f}'.format(bias))
+    return bias
 
 
 # Command line parsing
@@ -33,6 +31,10 @@ if options.list_cases:
     exit()
 
 # Run specific case, if requested
+case_nums = []
+case_names = []
+biases = []
+codes = []
 if options.case_name:
     if options.case_name in repo.names:
         case_num = repo.names[options.case_name]
@@ -44,10 +46,15 @@ if options.case_name:
         print('Running Case\n')
         case.make_model()
         code = case.execute(False)
-        if code != 0:
-            print(case.name + ' Failed Execution!')
+        if code:
+            bias = print_case(case)
         else:
-            print_case(case)
+            print(case.name + ' Failed Execution!')
+            bias = 0.
+        biases.append(bias)
+        case_names.append(case.name)
+        case_nums.append(case_num)
+        codes.append(code)
 
     else:
         print('Invalid Case Name: ' + options.case)
@@ -62,7 +69,26 @@ else:
         print(case.number, case.name)
         case.make_model()
         code = case.execute(True)
-        if code != 0:
-            print(case.name + ' Failed Execution!')
+        if code:
+            bias = print_case(case)
         else:
-            print_case(case)
+            print(case.name + ' Failed Execution!')
+            bias = 0.
+        biases.append(bias)
+        case_names.append(case.name)
+        codes.append(code)
+        case_nums.append(case.number)
+
+# Write to a CSV
+with open("results.csv", mode="w") as csv_file:
+    writer = csv.writer(csv_file, delimiter=",", quotechar='"',
+                        quoting=csv.QUOTE_NONNUMERIC)
+
+    # write the header
+    writer.writerow(['Case ID', 'Case Name', 'Eigenvalue Bias [pcm]',
+                     'Eigenvalue Bias Std. Dev. [pcm]', 'Successful Execution'])
+
+    # Now run each case
+    for i in range(len(case_nums)):
+        writer.writerow([case_nums[i], case_names[i], biases[i].n,
+                         biases[i].s, codes[i]])