diff --git a/wetb/prepost/h2_vs_hs2.py b/wetb/prepost/h2_vs_hs2.py
index b9d415786530c0118b512fc3c7d3aecf4ffe57d2..cd065486d076d1bd168c9183202cd2e47eaee165 100644
--- a/wetb/prepost/h2_vs_hs2.py
+++ b/wetb/prepost/h2_vs_hs2.py
@@ -442,7 +442,7 @@ class Sims(object):
         """
         Read a HAWCStab2 controller tuning file and return as tags
         """
-        tuning = hs2.hs2_control_tuning()
+        tuning = hs2.ReadControlTuning()
         tuning.read_parameters(fpath)
 
         tune_tags = {}
diff --git a/wetb/prepost/hawcstab2.py b/wetb/prepost/hawcstab2.py
index 00f705e93d8c53389d4860378ad0971ca40b2b54..1129dc7e0daf1299454975527580ef9ac5e0f3be 100644
--- a/wetb/prepost/hawcstab2.py
+++ b/wetb/prepost/hawcstab2.py
@@ -15,9 +15,6 @@ from future import standard_library
 standard_library.install_aliases()
 from builtins import object
 
-
-
-import unittest
 import os
 import re
 
@@ -28,8 +25,8 @@ from wetb.prepost import mplutils
 
 
 class dummy(object):
-    def __init__(self):
-        pass
+    def __init__(self, name='dummy'):
+        self.__name__ = name
 
 
 def ReadFileHAWCStab2Header(fname, widths):
@@ -172,6 +169,66 @@ class results(object):
                     'T_aero']
             self.operation = pd.DataFrame(operation, columns=cols)
 
+    def load_matrices(self, fpath, basename, operating_point=1,
+                      control_mat=False, local_wind_mat=False):
+        """Load HAWCStab2 State Space system matrices
+
+        The general file name format is:
+        BASENAMETYPE_ase_ops_OPERATING_POINT_NUMBER.dat
+
+        Where TYPE can be of the following:
+            * amat, bmat, bvmat, cmat, dmat, dvmat, emat, fmat, fvmat
+
+        Additionally, when including the control matrices:
+            * BASENAMETYPE_ops_OPERATING_POINT_NUMBER.dat
+            * TYPE: acmat, bcmat, ccmat, dcmat
+
+        Or when including local wind speed
+            * BASENAMETYPE_ase_ops_OPERATING_POINT_NUMBER.dat
+            * TYPE: bvmat_loc_v, dvmat_loc_v, fvmat_loc_v
+
+
+        Parameters
+        ----------
+
+        fpath : str
+
+        basename : str
+
+        operating_point : int, default=1
+
+
+        Returns
+        -------
+
+        matrices : dict
+
+        """
+        mnames = ['amat', 'bmat', 'bvmat', 'cmat', 'dmat', 'dvmat', 'emat',
+                  'fmat', 'fvmat']
+        mnames_c = ['acmat', 'bcmat', 'ccmat', 'dcmat']
+        mnames_v = ['bvmat_loc_v', 'dvmat_loc_v', 'fvmat_loc_v']
+
+        if control_mat:
+            mnames += mnames_c
+        if local_wind_mat:
+            mnames += mnames_v
+
+        matrices = {}
+
+        ase_template = '{:s}{:s}_ase_ops_{:d}.dat'
+        ops_template = '{:s}{:s}_ops_{:d}.dat'
+
+        for mname in mnames:
+            rpl = (basename, mname, operating_point)
+            template = ase_template
+            if mname in mnames_c:
+                template = ops_template
+            fname = os.path.join(fpath, template.format(*rpl))
+            matrices[mname] = np.loadtxt(fname)
+
+        return matrices
+
     def write_ae_sections_h2(self):
         """
         Get the aerosection positions from the HS2 ind result file and
@@ -237,7 +294,7 @@ class results(object):
         print('done!')
 
 
-class hs2_control_tuning(object):
+class ReadControlTuning(object):
 
     def __init__(self):
         """
@@ -293,4 +350,4 @@ class hs2_control_tuning(object):
 
 if __name__ == '__main__':
 
-    unittest.main()
+    pass
diff --git a/wetb/prepost/tests/test_hawcstab2.py b/wetb/prepost/tests/test_hawcstab2.py
index 249df2097c7a4167640d420eb3af1a1064981151..90d623ef6df16f587eebc33ba7ad89bb47cca78f 100644
--- a/wetb/prepost/tests/test_hawcstab2.py
+++ b/wetb/prepost/tests/test_hawcstab2.py
@@ -15,7 +15,7 @@ import os
 
 import numpy as np
 
-from wetb.prepost.hawcstab2 import results, hs2_control_tuning
+from wetb.prepost.hawcstab2 import results, ReadControlTuning
 
 
 class Tests(unittest.TestCase):
@@ -49,7 +49,7 @@ class Tests(unittest.TestCase):
 
     def test_linear_file(self):
 
-        hs2 = hs2_control_tuning()
+        hs2 = ReadControlTuning()
         hs2.read_parameters(self.fpath_linear)
 
         self.assertEqual(hs2.pi_gen_reg1.K, 0.108313E+07)
@@ -69,7 +69,7 @@ class Tests(unittest.TestCase):
 
     def test_quadratic_file(self):
 
-        hs2 = hs2_control_tuning()
+        hs2 = ReadControlTuning()
         hs2.read_parameters(self.fpath_quad)
 
         self.assertEqual(hs2.pi_gen_reg1.K, 0.108313E+07)