diff --git a/example/sofast_fringe/example_view_camera_distortion.py b/example/camera_calibration/example_view_camera_distortion.py
similarity index 100%
rename from example/sofast_fringe/example_view_camera_distortion.py
rename to example/camera_calibration/example_view_camera_distortion.py
diff --git a/example/sofast_fringe/example_multi_facet_data_process.py b/example/sofast_fringe/example_process_facet_ensemble.py
similarity index 83%
rename from example/sofast_fringe/example_multi_facet_data_process.py
rename to example/sofast_fringe/example_process_facet_ensemble.py
index 0f2f78e85..b536549aa 100644
--- a/example/sofast_fringe/example_multi_facet_data_process.py
+++ b/example/sofast_fringe/example_process_facet_ensemble.py
@@ -1,25 +1,27 @@
-import os
-from os.path import join
+from os.path import join, dirname
 
 import matplotlib
 
+from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
+from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble
+from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
 from opencsp.app.sofast.lib.ImageCalibrationScaling import ImageCalibrationScaling
 from opencsp.app.sofast.lib.MeasurementSofastFringe import (
     MeasurementSofastFringe as Measurement,
 )
 from opencsp.app.sofast.lib.ProcessSofastFringe import ProcessSofastFringe as Sofast
+from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
 from opencsp.app.sofast.lib.visualize_setup import visualize_setup
 from opencsp.common.lib.camera.Camera import Camera
 from opencsp.common.lib.csp.FacetEnsemble import FacetEnsemble
-from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
-from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble
-from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
-from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
+from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
 import opencsp.common.lib.render.figure_management as fm
 import opencsp.common.lib.render_control.RenderControlAxis as rca
 import opencsp.common.lib.render_control.RenderControlFigure as rcfg
 import opencsp.common.lib.render_control.RenderControlMirror as rcm
+import opencsp.common.lib.tool.log_tools as lt
+import opencsp.common.lib.tool.file_tools as ft
 
 
 def example_driver():
@@ -44,9 +46,11 @@ def example_driver():
     file_ensemble = join(sample_data_dir, 'Ensemble_lab_6x4.json')
 
     # Define save dir
-    dir_save = join(os.path.dirname(__file__), 'data/output/facet_ensemble')
-    if not os.path.exists(dir_save):
-        os.makedirs(dir_save)
+    dir_save = join(dirname(__file__), 'data/output/facet_ensemble')
+    ft.create_directories_if_necessary(dir_save)
+
+    # Set up logger
+    lt.logger(join(dir_save, 'log.txt'))
 
     # Load data
     camera = Camera.load_from_hdf(file_camera)
@@ -58,12 +62,9 @@ def example_driver():
     # Define facet data
     facet_data = [DefinitionFacet.load_from_json(file_facet)] * ensemble_data.num_facets
 
-    # Define surface data (plano)
-    # surface_data = [dict(surface_type='plano', robust_least_squares=False, downsample=20)] * ensemble_data.num_facets
-    # Define surface data (parabolic)
+    # Define surface data
     surface_data = [
-        dict(
-            surface_type='parabolic',
+        Surface2DParabolic(
             initial_focal_lengths_xy=(100.0, 100.0),
             robust_least_squares=False,
             downsample=20,
@@ -88,12 +89,10 @@ def example_driver():
 
     # Calculate focal length from parabolic fit
     for idx in range(sofast.num_facets):
-        if surface_data[idx]['surface_type'] == 'parabolic':
-            surf_coefs = sofast.data_characterization_facet[idx].surf_coefs_facet
-            focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]]
-            print('Parabolic fit focal lengths:')
-            print(f'  X {focal_lengths_xy[0]:.3f} m')
-            print(f'  Y {focal_lengths_xy[1]:.3f} m')
+        surf_coefs = sofast.data_characterization_facet[idx].surf_coefs_facet
+        focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]]
+        lt.info(f'Facet {idx:d} xy focal lengths (meters): '
+                f'{focal_lengths_xy[0]:.3f}, {focal_lengths_xy[1]:.3f}')
 
     # Get optic representation
     ensemble: FacetEnsemble = sofast.get_optic()
diff --git a/example/sofast_fringe/example_single_facet_data_process.py b/example/sofast_fringe/example_process_single_facet.py
similarity index 79%
rename from example/sofast_fringe/example_single_facet_data_process.py
rename to example/sofast_fringe/example_process_single_facet.py
index fc062b9a1..393a526d6 100644
--- a/example/sofast_fringe/example_single_facet_data_process.py
+++ b/example/sofast_fringe/example_process_single_facet.py
@@ -1,23 +1,25 @@
-import os
-from os.path import join
+from os.path import join, dirname
 
 import matplotlib
 
+from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
+from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
 from opencsp.app.sofast.lib.ImageCalibrationScaling import ImageCalibrationScaling
 from opencsp.app.sofast.lib.MeasurementSofastFringe import (
     MeasurementSofastFringe as Measurement,
 )
 from opencsp.app.sofast.lib.ProcessSofastFringe import ProcessSofastFringe as Sofast
+from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
 from opencsp.app.sofast.lib.visualize_setup import visualize_setup
 from opencsp.common.lib.camera.Camera import Camera
 from opencsp.common.lib.csp.Facet import Facet
-from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
-from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
-from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
+from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
 import opencsp.common.lib.render.figure_management as fm
 import opencsp.common.lib.render_control.RenderControlAxis as rca
 import opencsp.common.lib.render_control.RenderControlFigure as rcfg
+import opencsp.common.lib.tool.log_tools as lt
+import opencsp.common.lib.tool.file_tools as ft
 
 
 def example_driver():
@@ -40,9 +42,11 @@ def example_driver():
     file_facet = join(sample_data_dir, 'Facet_NSTTF.json')
 
     # Define save dir
-    dir_save = join(os.path.dirname(__file__), 'data/output/single_facet')
-    if not os.path.exists(dir_save):
-        os.makedirs(dir_save)
+    dir_save = join(dirname(__file__), 'data/output/single_facet')
+    ft.create_directories_if_necessary(dir_save)
+
+    # Set up logger
+    lt.logger(join(dir_save, 'log.txt'))
 
     # Load data
     camera = Camera.load_from_hdf(file_camera)
@@ -52,9 +56,8 @@ def example_driver():
     facet_data = DefinitionFacet.load_from_json(file_facet)
 
     # Define surface definition (parabolic surface)
-    surface_data = dict(
-        surface_type='parabolic',
-        initial_focal_lengths_xy=(300.0, 300),
+    surface = Surface2DParabolic(
+        initial_focal_lengths_xy=(300., 300.),
         robust_least_squares=True,
         downsample=10,
     )
@@ -66,15 +69,13 @@ def example_driver():
     sofast = Sofast(measurement, camera, display)
 
     # Process
-    sofast.process_optic_singlefacet(facet_data, surface_data)
+    sofast.process_optic_singlefacet(facet_data, surface)
 
     # Calculate focal length from parabolic fit
-    if surface_data['surface_type'] == 'parabolic':
-        surf_coefs = sofast.data_characterization_facet[0].surf_coefs_facet
-        focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]]
-        print('Parabolic fit focal lengths:')
-        print(f'  X {focal_lengths_xy[0]:.3f} m')
-        print(f'  Y {focal_lengths_xy[1]:.3f} m')
+    surf_coefs = sofast.data_characterization_facet[0].surf_coefs_facet
+    focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]]
+    lt.info(f'Facet xy focal lengths (meters): '
+            f'{focal_lengths_xy[0]:.3f}, {focal_lengths_xy[1]:.3f}')
 
     # Get optic representation
     facet: Facet = sofast.get_optic()
diff --git a/example/sofast_fringe/example_undefined_facet_data_process.py b/example/sofast_fringe/example_process_undefined_shape.py
similarity index 78%
rename from example/sofast_fringe/example_undefined_facet_data_process.py
rename to example/sofast_fringe/example_process_undefined_shape.py
index cb8c71f90..212d69e9a 100644
--- a/example/sofast_fringe/example_undefined_facet_data_process.py
+++ b/example/sofast_fringe/example_process_undefined_shape.py
@@ -1,5 +1,6 @@
-import os
-from os.path import join
+from os.path import join, dirname
+
+import matplotlib
 
 from opencsp.app.sofast.lib.visualize_setup import visualize_setup
 from opencsp.app.sofast.lib.ImageCalibrationScaling import ImageCalibrationScaling
@@ -11,10 +12,13 @@
 from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
 from opencsp.common.lib.camera.Camera import Camera
 from opencsp.common.lib.csp.Facet import Facet
+from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
 import opencsp.common.lib.render.figure_management as fm
 import opencsp.common.lib.render_control.RenderControlAxis as rca
 import opencsp.common.lib.render_control.RenderControlFigure as rcfg
+import opencsp.common.lib.tool.log_tools as lt
+import opencsp.common.lib.tool.file_tools as ft
 
 
 def example_driver():
@@ -36,9 +40,11 @@ def example_driver():
     file_calibration = join(sample_data_dir, 'image_calibration.h5')
 
     # Save directory
-    dir_save = join(os.path.dirname(__file__), 'data/output/undefined_facet')
-    if not os.path.exists(dir_save):
-        os.makedirs(dir_save)
+    dir_save = join(dirname(__file__), 'data/output/undefined_facet')
+    ft.create_directories_if_necessary(dir_save)
+
+    # Set up logger
+    lt.logger(join(dir_save, 'log.txt'))
 
     # Load data
     camera = Camera.load_from_hdf(file_camera)
@@ -47,9 +53,8 @@ def example_driver():
     calibration = ImageCalibrationScaling.load_from_hdf(file_calibration)
 
     # Define surface definition (parabolic surface)
-    surface_data = dict(
-        surface_type='parabolic',
-        initial_focal_lengths_xy=(300.0, 300),
+    surface = Surface2DParabolic(
+        initial_focal_lengths_xy=(300., 300.),
         robust_least_squares=True,
         downsample=10,
     )
@@ -62,15 +67,13 @@ def example_driver():
     sofast.params.mask_keep_largest_area = True
 
     # Process
-    sofast.process_optic_undefined(surface_data)
+    sofast.process_optic_undefined(surface)
 
     # Calculate focal length from parabolic fit
-    if surface_data['surface_type'] == 'parabolic':
-        surf_coefs = sofast.data_characterization_facet[0].surf_coefs_facet
-        focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]]
-        print('Parabolic fit focal lengths:')
-        print(f'  X {focal_lengths_xy[0]:.3f} m')
-        print(f'  Y {focal_lengths_xy[1]:.3f} m')
+    surf_coefs = sofast.data_characterization_facet[0].surf_coefs_facet
+    focal_lengths_xy = [1 / 4 / surf_coefs[2], 1 / 4 / surf_coefs[5]]
+    lt.info(f'Facet xy focal lengths (meters): '
+            f'{focal_lengths_xy[0]:.3f}, {focal_lengths_xy[1]:.3f}')
 
     # Get optic representation
     facet: Facet = sofast.get_optic()
diff --git a/opencsp/app/sofast/calibration/test/data/output/CalibrationScreenPosition/screen_distortion_data_100_100.h5 b/opencsp/app/sofast/calibration/test/data/output/CalibrationScreenPosition/screen_distortion_data_100_100.h5
deleted file mode 100644
index 45bd931f6..000000000
Binary files a/opencsp/app/sofast/calibration/test/data/output/CalibrationScreenPosition/screen_distortion_data_100_100.h5 and /dev/null differ
diff --git a/opencsp/app/sofast/lib/ProcessSofastFringe.py b/opencsp/app/sofast/lib/ProcessSofastFringe.py
index c8790bd0f..bbace7091 100644
--- a/opencsp/app/sofast/lib/ProcessSofastFringe.py
+++ b/opencsp/app/sofast/lib/ProcessSofastFringe.py
@@ -3,36 +3,35 @@
 """
 
 from typing import Literal
-import warnings
 
 import numpy as np
 
-from opencsp.app.sofast.lib.MeasurementSofastFringe import (
-    MeasurementSofastFringe as Measurement,
-)
-from opencsp.app.sofast.lib.ParamsSofastFringe import ParamsSofastFringe
-from opencsp.common.lib.camera.Camera import Camera
-from opencsp.common.lib.csp.Facet import Facet
-from opencsp.common.lib.csp.FacetEnsemble import FacetEnsemble
-from opencsp.common.lib.csp.MirrorPoint import MirrorPoint
 import opencsp.app.sofast.lib.calculation_data_classes as cdc
-from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
 from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble
 from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
+from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
 import opencsp.app.sofast.lib.image_processing as ip
+from opencsp.app.sofast.lib.MeasurementSofastFringe import MeasurementSofastFringe as Measurement
+from opencsp.app.sofast.lib.ParamsSofastFringe import ParamsSofastFringe
 import opencsp.app.sofast.lib.process_optics_geometry as po
+from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
+from opencsp.common.lib.camera.Camera import Camera
+from opencsp.common.lib.csp.Facet import Facet
+from opencsp.common.lib.csp.FacetEnsemble import FacetEnsemble
+from opencsp.common.lib.csp.MirrorPoint import MirrorPoint
 from opencsp.common.lib.deflectometry.SlopeSolver import SlopeSolver
 from opencsp.common.lib.deflectometry.SlopeSolverData import SlopeSolverData
-from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
+from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract
 from opencsp.common.lib.geometry.RegionXY import RegionXY
 from opencsp.common.lib.geometry.TransformXYZ import TransformXYZ
 from opencsp.common.lib.geometry.Uxyz import Uxyz
 from opencsp.common.lib.geometry.Vxy import Vxy
 from opencsp.common.lib.geometry.Vxyz import Vxyz
-from opencsp.common.lib.tool.hdf5_tools import save_hdf5_datasets
+from opencsp.common.lib.tool.hdf5_tools import HDF5_SaveAbstract
+import opencsp.common.lib.tool.log_tools as lt
 
 
-class ProcessSofastFringe:
+class ProcessSofastFringe(HDF5_SaveAbstract):
     """Class that processes measurement data captured by a SOFAST
     system. Computes optic surface slope and saves data to HDF5 format.
 
@@ -192,7 +191,7 @@ def __init__(
         self.data_facet_def: list[DefinitionFacet] = None
         self.data_ensemble_def: DefinitionEnsemble = None
 
-        self.data_surface_params: list[dict] = None
+        self.data_surfaces: list[Surface2DAbstract] = None
 
         self.data_geometry_general: cdc.CalculationDataGeometryGeneral = None
         self.data_image_processing_general: cdc.CalculationImageProcessingGeneral = None
@@ -209,52 +208,16 @@ def help(self) -> None:
         """Prints Sofast doc string"""
         print(self.__doc__)
 
-    @staticmethod
-    def _check_surface_data(surf_data: dict) -> None:
-        """Checks that all necessary fields are present in surface data dict"""
-        if 'surface_type' not in surf_data.keys():
-            raise ValueError('Missing "surface_type" key in surface_data dictionary.')
-
-        if surf_data['surface_type'] == 'parabolic':
-            fields_exp = [
-                'surface_type',
-                'initial_focal_lengths_xy',
-                'robust_least_squares',
-                'downsample',
-            ]
-        elif surf_data['surface_type'] == 'plano':
-            fields_exp = ['surface_type', 'robust_least_squares', 'downsample']
-        else:
-            raise ValueError(
-                f'Given surface type {surf_data["surface_type"]} is not supported.'
-            )
-
-        for k in surf_data.keys():
-            if k in fields_exp:
-                idx = fields_exp.index(k)
-                fields_exp.pop(idx)
-            else:
-                raise ValueError(
-                    f'Unrecognized field, {k}, in surface_data dictionary.'
-                )
-
-        if len(fields_exp) > 0:
-            raise ValueError(f'Missing fields in surface_data dictionary: {fields_exp}')
-
-    def process_optic_undefined(self, surface_data: dict) -> None:
+    def process_optic_undefined(self, surface: Surface2DAbstract) -> None:
         """
         Processes optic geometry, screen intersection points, and solves
         for slops for undefined optical surface.
 
         Parameters
         ----------
-        surface_data : dict
-            See Sofast documentation or Sofast.help() for more details.
-
+        surface_data : Surface2DAbstract
+            Surface type definition
         """
-        # Check input data
-        self._check_surface_data(surface_data)
-
         # Process optic/setup geometry
         self._process_optic_undefined_geometry()
 
@@ -262,11 +225,11 @@ def process_optic_undefined(self, surface_data: dict) -> None:
         self._process_display()
 
         # Solve slopes
-        self._solve_slopes([surface_data])
+        self._solve_slopes([surface])
 
-    def process_optic_singlefacet(
-        self, facet_data: DefinitionFacet, surface_data: dict
-    ) -> None:
+    def process_optic_singlefacet(self,
+                                  facet_data: DefinitionFacet,
+                                  surface: Surface2DAbstract) -> None:
         """
         Processes optic geometry, screen intersection points, and solves
         for slops for single facet optic.
@@ -275,13 +238,9 @@ def process_optic_singlefacet(
         ----------
         facet_data : DefinitionFacet
             Facet data object.
-        surface_data : dict
-            See Sofast documentation or Sofast.help() for more details.
-
+        surface_data : Surface2DAbstract
+            Surface type definition.
         """
-        # Check input data
-        self._check_surface_data(surface_data)
-
         # Process optic/setup geometry
         self._process_optic_singlefacet_geometry(facet_data)
 
@@ -289,14 +248,12 @@ def process_optic_singlefacet(
         self._process_display()
 
         # Solve slopes
-        self._solve_slopes([surface_data])
+        self._solve_slopes([surface])
 
-    def process_optic_multifacet(
-        self,
-        facet_data: list[DefinitionFacet],
-        ensemble_data: DefinitionEnsemble,
-        surface_data: list[dict],
-    ) -> None:
+    def process_optic_multifacet(self,
+                                 facet_data: list[DefinitionFacet],
+                                 ensemble_data: DefinitionEnsemble,
+                                 surfaces: list[Surface2DAbstract]) -> None:
         """
         Processes optic geometry, screen intersection points, and solves
         for slops for multi-facet optic.
@@ -312,11 +269,11 @@ def process_optic_multifacet(
 
         """
         # Check inputs
-        if len(facet_data) != len(surface_data):
+        if len(facet_data) != len(surfaces):
             raise ValueError(
-                f'Length of facet_data does not equal length of surface data; facet_data={len(facet_data)}, surface_data={len(surface_data)}'
+                'Length of facet_data does not equal length of surfaces'
+                f'facet_data={len(facet_data)}, surface_data={len(surfaces)}'
             )
-        list(map(self._check_surface_data, surface_data))
 
         # Process optic/setup geometry
         self._process_optic_multifacet_geometry(facet_data, ensemble_data)
@@ -325,7 +282,7 @@ def process_optic_multifacet(
         self._process_display()
 
         # Solve slopes
-        self._solve_slopes(surface_data)
+        self._solve_slopes(surfaces)
 
         # Calculate facet pointing
         self._calculate_facet_pointing()
@@ -383,9 +340,9 @@ def _process_optic_singlefacet_geometry(self, facet_data: DefinitionFacet) -> No
         self.optic_type = 'single'
 
         if self.params.geometry_data_debug.debug_active:
-            print('Sofast image processing debug on.')
+            lt.info('Sofast image processing debug on.')
         if self.params.slope_solver_data_debug.debug_active:
-            print('SlopeSolver debug on.')
+            lt.info('SlopeSolver debug on.')
 
         # Calculate raw mask
         params = [
@@ -425,8 +382,7 @@ def _process_optic_singlefacet_geometry(self, facet_data: DefinitionFacet) -> No
     def _process_optic_multifacet_geometry(
         self, facet_data: list[DefinitionFacet], ensemble_data: DefinitionEnsemble
     ) -> None:
-        """
-        Processes optic geometry for an ensemble of facets.
+        """Processes optic geometry for an ensemble of facets.
 
         Parameters
         ----------
@@ -434,7 +390,6 @@ def _process_optic_multifacet_geometry(
             List of DefinitionFacet objects.
         ensemble_data : DefinitionEnsemble
             Ensemble data object.
-
         """
         # Get number of facets
         self.num_facets = ensemble_data.num_facets
@@ -443,7 +398,8 @@ def _process_optic_multifacet_geometry(
         # Check inputs
         if len(facet_data) != self.num_facets:
             raise ValueError(
-                f'Given length of facet data is {len(facet_data):d} but ensemble_data expects {ensemble_data.num_facets:d} facets.'
+                f'Given length of facet data is {len(facet_data):d}'
+                f'but ensemble_data expects {ensemble_data.num_facets:d} facets.'
             )
 
         # Calculate mask
@@ -456,9 +412,9 @@ def _process_optic_multifacet_geometry(
         mask_raw = ip.calc_mask_raw(self.measurement.mask_images, *params)
 
         if self.params.mask_keep_largest_area:
-            warnings.warn(
-                '"keep_largest_area" mask processing option cannot be used for multifacet ensembles. This will be turned off.',
-                stacklevel=2,
+            lt.warn(
+                '"keep_largest_area" mask processing option cannot be used '
+                'for multifacet ensembles. This will be turned off.',
             )
             self.params.mask_keep_largest_area = False
 
@@ -523,9 +479,9 @@ def _process_display(self) -> None:
             nan_mask = np.isnan(v_screen_points_screen.data).sum(0).astype(bool)
             mask_bad_pixels = np.zeros(mask_processed.shape, dtype=bool)
             if np.any(nan_mask):
-                warnings.warn(
-                    f'{nan_mask.sum():d} / {nan_mask.size:d} points are NANs in calculated screen points for facet {idx_facet:d}. These data points will be removed.',
-                    stacklevel=2,
+                lt.warn(
+                    f'{nan_mask.sum():d} / {nan_mask.size:d} points are NANs in calculated '
+                    'screen points for facet {idx_facet:d}. These data points will be removed.',
                 )
                 # Make mask of NANs
                 mask_bad_pixels[mask_processed] = nan_mask
@@ -563,16 +519,14 @@ def _process_display(self) -> None:
                 v_screen_points_facet
             )
 
-    def _solve_slopes(self, surface_data: list[dict]) -> None:
+    def _solve_slopes(self, surfaces: list[Surface2DAbstract]) -> None:
         """
         Solves slopes of each active pixel for each facet.
 
         Parameters
         ----------
-        surface_data : list[dict]
-            List containing one dictionary for each facet being processed.
-            See SlopeSolver documentation for details.
-
+        surface_data : list[Surface2DAbstract]
+            List of surface definition classes.
         """
         # Check inputs
         if self.data_geometry_facet is None:
@@ -610,7 +564,7 @@ def _solve_slopes(self, surface_data: list[dict]) -> None:
                 'dist_optic_screen': self.data_geometry_facet[
                     facet_idx
                 ].measure_point_screen_distance,
-                'surface_data': surface_data[facet_idx],
+                'surface': surfaces[facet_idx],
                 'debug': self.params.slope_solver_data_debug,
             }
 
@@ -627,7 +581,7 @@ def _solve_slopes(self, surface_data: list[dict]) -> None:
             self.data_characterization_facet.append(slope_solver.get_data())
 
         # Save input surface parameters data
-        self.data_surface_params = [s.copy() for s in surface_data]
+        self.data_surfaces = surfaces
 
     def _calculate_facet_pointing(
         self, reference: Literal['average'] | int = 'average'
@@ -785,69 +739,58 @@ def get_optic(
         else:
             return facets[0]
 
-    def save_to_hdf(self, file: str) -> None:
-        """
-        Saves all processed data to HDF file. See class docstring
-        for more information on data output format.
+    def save_to_hdf(self, file: str, prefix: str = ''):
+        """Saves data to given file. Data is stored as: PREFIX + Folder/Field_1
 
         Parameters
         ----------
         file : str
-            HDF file name.
+            HDF file to save to
+        prefix : str
+            Prefix to append to folder path within HDF file (folders must be separated by "/")
         """
+        # Log
+        lt.info(f'Saving SofastFringe data to: {file:s}, in HDF5 folder: "{prefix:s}"')
+
         # One per measurement
         if self.data_error is not None:
-            self.data_error.save_to_hdf(file, 'DataSofastCalculation/general/')
-        self.data_geometry_general.save_to_hdf(file, 'DataSofastCalculation/general/')
+            self.data_error.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/')
+        self.data_geometry_general.save_to_hdf(file, f'{prefix:s}DataSofastCalculation/general/')
         self.data_image_processing_general.save_to_hdf(
-            file, 'DataSofastCalculation/general/'
-        )
+            file, f'{prefix:s}DataSofastCalculation/general/')
 
         # Sofast parameters
-        self.params.save_to_hdf(file, 'DataSofastInput/')
+        self.params.save_to_hdf(file, f'{prefix:s}DataSofastInput/')
 
         # Facet definition
         if self.data_facet_def is not None:
             for idx_facet, facet_data in enumerate(self.data_facet_def):
                 facet_data.save_to_hdf(
-                    file, f'DataSofastInput/optic_definition/facet_{idx_facet:03d}/'
-                )
+                    file, f'{prefix:s}DataSofastInput/optic_definition/facet_{idx_facet:03d}/')
 
         # Ensemble definition
         if self.data_ensemble_def is not None:
-            self.data_ensemble_def.save_to_hdf(
-                file, 'DataSofastInput/optic_definition/'
-            )
+            self.data_ensemble_def.save_to_hdf(file, f'{prefix:s}DataSofastInput/optic_definition/')
 
         # Surface definition
-        # TODO: make surface_params a data class
-        for idx_facet, surface_params in enumerate(self.data_surface_params):
-            data = list(surface_params.values())
-            datasets = list(surface_params.keys())
-            datasets = [
-                f'DataSofastInput/optic_definition/facet_{idx_facet:03d}/surface_definition/'
-                + d
-                for d in datasets
-            ]
-            save_hdf5_datasets(data, datasets, file)
+        for idx_facet, surface in enumerate(self.data_surfaces):
+            surface.save_to_hdf(
+                file,
+                f'{prefix:s}DataSofastInput/optic_definition/facet_{idx_facet:03d}/')
 
         # Calculations, one per facet
         for idx_facet in range(self.num_facets):
             # Save facet slope data
             self.data_characterization_facet[idx_facet].save_to_hdf(
-                file, f'DataSofastCalculation/facet/facet_{idx_facet:03d}/'
-            )
+                file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/')
             # Save facet geometry data
             self.data_geometry_facet[idx_facet].save_to_hdf(
-                file, f'DataSofastCalculation/facet/facet_{idx_facet:03d}/'
-            )
+                file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/')
             # Save facet image processing data
             self.data_image_processing_facet[idx_facet].save_to_hdf(
-                file, f'DataSofastCalculation/facet/facet_{idx_facet:03d}/'
-            )
+                file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/')
 
             if self.data_characterization_ensemble:
                 # Save ensemle data
                 self.data_characterization_ensemble[idx_facet].save_to_hdf(
-                    file, f'DataSofastCalculation/facet/facet_{idx_facet:03d}/'
-                )
+                    file, f'{prefix:s}DataSofastCalculation/facet/facet_{idx_facet:03d}/')
diff --git a/opencsp/app/sofast/lib/process_optics_geometry.py b/opencsp/app/sofast/lib/process_optics_geometry.py
index ed2597954..efa98adf3 100644
--- a/opencsp/app/sofast/lib/process_optics_geometry.py
+++ b/opencsp/app/sofast/lib/process_optics_geometry.py
@@ -61,11 +61,18 @@ def process_singlefacet_geometry(
 
     Returns
     -------
-    calculation_data_classes.CalculationDataGeometryGeneral
-    calculation_data_classes.CalculationImageProcessingGeneral
-    list[calculation_data_classes.CalculationDataGeometryFacet]
-    list[calculation_data_classes.CalculationImageProcessingFacet]
-    calculation_data_classes.CalculationError]
+    data_geometry_general: calculation_data_classes.CalculationDataGeometryGeneral
+        Positional optic geometry calculations general to entire measurement; not facet specific.
+    data_image_processing_general: calculation_data_classes.CalculationImageProcessingGeneral
+        Image processing calculations general to entire measurement; not facet specific.
+    data_geometry_facet: list[calculation_data_classes.CalculationDataGeometryFacet]
+        List of positional optic geometry calculations specific to each facet. Order is
+        same as input facet definitions.
+    data_image_processing_facet: list[calculation_data_classes.CalculationImageProcessingFacet]
+        List of image processing calcualtions specific to each facet. Order is same as input facet
+        definitions.
+    data_error: calculation_data_classes.CalculationError
+        Geometric/positional errors and reprojection errors associated with solving for facet location.
     """
     if debug.debug_active:
         print('Image processing debug on.')
@@ -305,7 +312,7 @@ def process_undefined_geometry(
     list[cdc.CalculationImageProcessingFacet],
     cdc.CalculationError,
 ]:
-    """Processes optic geometry for undefined deflectrometry measurement
+    """Processes optic geometry for undefined deflectometry measurement
 
     Parameters
     ----------
@@ -324,11 +331,18 @@ def process_undefined_geometry(
 
     Returns
     -------
-    calculation_data_classes.CalculationDataGeometryGeneral
-    calculation_data_classes.CalculationImageProcessingGeneral
-    list[calculation_data_classes.CalculationDataGeometryFacet]
-    list[calculation_data_classes.CalculationImageProcessingFacet]
-    calculation_data_classes.CalculationError]
+    data_geometry_general: calculation_data_classes.CalculationDataGeometryGeneral
+        Positional optic geometry calculations general to entire measurement; not facet specific.
+    data_image_processing_general: calculation_data_classes.CalculationImageProcessingGeneral
+        Image processing calculations general to entire measurement; not facet specific.
+    data_geometry_facet: list[calculation_data_classes.CalculationDataGeometryFacet]
+        List of positional optic geometry calculations specific to each facet. Order is
+        same as input facet definitions.
+    data_image_processing_facet: list[calculation_data_classes.CalculationImageProcessingFacet]
+        List of image processing calcualtions specific to each facet. Order is same as input facet
+        definitions.
+    data_error: calculation_data_classes.CalculationError
+        Geometric/positional errors and reprojection errors associated with solving for facet location.
     """
     if debug.debug_active:
         print(
@@ -409,6 +423,44 @@ def process_multifacet_geometry(
     list[cdc.CalculationImageProcessingFacet],
     cdc.CalculationError,
 ]:
+    """Processes optic geometry for multifacet deflectometry measurement
+
+    Parameters
+    ----------
+    facet_data : DefinitionFacet
+        Facet definition object
+    ensemble_data : DefinitionEnsemble
+        Ensemble definition object
+    mask_raw : ndarray
+        Raw calculated mask, shape (m, n) array of booleans
+    v_meas_pt_ensemble : Vxyz
+        Measure point lcoation on ensemble, meters
+    orientation : SpatialOrientation
+        SpatialOrientation object
+    camera : Camera
+        Camera object
+    optic_screen_dist : float
+        Optic to screen distance, meters
+    params : ParamsOpticGeometry, optional
+        ParamsOpticGeometry object, by default ParamsOpticGeometry()
+    debug : DebugOpticsGeometry, optional
+        DebugOpticsGeometry object, by default DebugOpticsGeometry()
+
+    Returns
+    -------
+    data_geometry_general: calculation_data_classes.CalculationDataGeometryGeneral
+        Positional optic geometry calculations general to entire measurement; not facet specific.
+    data_image_processing_general: calculation_data_classes.CalculationImageProcessingGeneral
+        Image processing calculations general to entire measurement; not facet specific.
+    data_geometry_facet: list[calculation_data_classes.CalculationDataGeometryFacet]
+        List of positional optic geometry calculations specific to each facet. Order is
+        same as input facet definitions.
+    data_image_processing_facet: list[calculation_data_classes.CalculationImageProcessingFacet]
+        List of image processing calcualtions specific to each facet. Order is same as input facet
+        definitions.
+    data_error: calculation_data_classes.CalculationError
+        Geometric/positional errors and reprojection errors associated with solving for facet location.
+    """
     if debug.debug_active:
         print(
             'Image processing debug on, but is not yet supported for undefined mirrors.'
diff --git a/opencsp/app/sofast/test/test_integration_multi_facet.py b/opencsp/app/sofast/test/test_integration_multi_facet.py
index 447c8cbef..bdf155592 100644
--- a/opencsp/app/sofast/test/test_integration_multi_facet.py
+++ b/opencsp/app/sofast/test/test_integration_multi_facet.py
@@ -7,15 +7,16 @@
 from scipy.spatial.transform import Rotation
 import numpy as np
 
+from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
+from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble
+from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
 from opencsp.app.sofast.lib.ImageCalibrationScaling import ImageCalibrationScaling
 from opencsp.app.sofast.lib.ProcessSofastFringe import ProcessSofastFringe as Sofast
 from opencsp.app.sofast.lib.MeasurementSofastFringe import (
     MeasurementSofastFringe as Measurement,
 )
 from opencsp.common.lib.camera.Camera import Camera
-from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
-from opencsp.app.sofast.lib.DefinitionEnsemble import DefinitionEnsemble
-from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
+from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
 from opencsp.common.lib.geometry.Vxyz import Vxyz
 from opencsp.common.lib.tool.hdf5_tools import load_hdf5_datasets
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
@@ -41,7 +42,6 @@ def setUpClass(cls, base_dir: str | None = None):
         # Directory Setup
         file_dataset = os.path.join(base_dir, 'calculations_facet_ensemble/data.h5')
         file_measurement = os.path.join(base_dir, 'measurement_ensemble.h5')
-        print(f'Using dataset: {os.path.abspath(file_dataset)}')
 
         # Load data
         camera = Camera.load_from_hdf(file_dataset)
@@ -123,20 +123,19 @@ def setUpClass(cls, base_dir: str | None = None):
             )
 
         # Load surface data
-        surface_data = []
+        surfaces = []
         for idx in range(len(facet_data)):
             datasets = [
                 f'DataSofastInput/surface_params/facet_{idx:03d}/downsample',
                 f'DataSofastInput/surface_params/facet_{idx:03d}/initial_focal_lengths_xy',
                 f'DataSofastInput/surface_params/facet_{idx:03d}/robust_least_squares',
-                f'DataSofastInput/surface_params/facet_{idx:03d}/surface_type',
             ]
             data = load_hdf5_datasets(datasets, file_dataset)
             data['robust_least_squares'] = bool(data['robust_least_squares'])
-            surface_data.append(data)
+            surfaces.append(Surface2DParabolic(**data))
 
         # Run SOFAST
-        sofast.process_optic_multifacet(facet_data, ensemble_data, surface_data)
+        sofast.process_optic_multifacet(facet_data, ensemble_data, surfaces)
 
         # Store data
         cls.data_test = {'slopes_facet_xy': [], 'surf_coefs_facet': []}
@@ -188,9 +187,4 @@ def test_surf_coefs(self):
 
 
 if __name__ == '__main__':
-    Test = TestMulti()
-    Test.setUpClass()
-
-    Test.test_slope()
-    Test.test_surf_coefs()
-    print('All tests run')
+    unittest.main()
diff --git a/opencsp/app/sofast/test/test_integration_single_facet.py b/opencsp/app/sofast/test/test_integration_single_facet.py
index 264688034..e5bf2a75c 100644
--- a/opencsp/app/sofast/test/test_integration_single_facet.py
+++ b/opencsp/app/sofast/test/test_integration_single_facet.py
@@ -7,14 +7,16 @@
 
 import numpy as np
 
+from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
+from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
 from opencsp.app.sofast.lib.ImageCalibrationScaling import ImageCalibrationScaling
 from opencsp.app.sofast.lib.MeasurementSofastFringe import (
     MeasurementSofastFringe as Measurement,
 )
 from opencsp.app.sofast.lib.ProcessSofastFringe import ProcessSofastFringe as Sofast
 from opencsp.common.lib.camera.Camera import Camera
-from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
-from opencsp.app.sofast.lib.DefinitionFacet import DefinitionFacet
+from opencsp.common.lib.deflectometry.Surface2DPlano import Surface2DPlano
+from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
 from opencsp.common.lib.geometry.Vxyz import Vxyz
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
 from opencsp.common.lib.tool.hdf5_tools import load_hdf5_datasets
@@ -43,10 +45,6 @@ def setUpClass(cls, base_dir: str | None = None):
         )
         if len(cls.files_dataset) == 0:
             raise ValueError('No single-facet datsets found.')
-        else:
-            print(f'Testing {len(cls.files_dataset)} single facet datasets')
-            for file in cls.files_dataset:
-                print(f'Using dataset: {os.path.abspath(file)}')
 
         # Define component files
         file_measurement = os.path.join(base_dir, 'measurement_facet.h5')
@@ -90,6 +88,16 @@ def setUpClass(cls, base_dir: str | None = None):
                         file_dataset,
                     )
                 )
+                surface = Surface2DParabolic(
+                    surface_data['initial_focal_lengths_xy'],
+                    surface_data['robust_least_squares'],
+                    surface_data['downsample'],
+                )
+            else:
+                surface = Surface2DPlano(
+                    surface_data['robust_least_squares'],
+                    surface_data['downsample'],
+                )
 
             # Load optic data
             facet_data = load_hdf5_datasets(
@@ -148,7 +156,7 @@ def setUpClass(cls, base_dir: str | None = None):
             ]
 
             # Run SOFAST
-            sofast.process_optic_singlefacet(facet_data, surface_data)
+            sofast.process_optic_singlefacet(facet_data, surface)
 
             # Store test data
             cls.slopes.append(sofast.data_characterization_facet[0].slopes_facet_xy)
@@ -191,16 +199,4 @@ def test_int_points(self):
 
 
 if __name__ == '__main__':
-    Test = TestSingle()
-    Test.setUpClass()
-
-    print('test_slopes', flush=True)
-    Test.test_slopes()
-
-    print('test_surf_coefs', flush=True)
-    Test.test_surf_coefs()
-
-    print('test_int_points', flush=True)
-    Test.test_int_points()
-
-    print('All tests run.')
+    unittest.main()
diff --git a/opencsp/app/sofast/test/test_integration_undefined.py b/opencsp/app/sofast/test/test_integration_undefined.py
index c0a3dfac0..856a67388 100644
--- a/opencsp/app/sofast/test/test_integration_undefined.py
+++ b/opencsp/app/sofast/test/test_integration_undefined.py
@@ -2,116 +2,109 @@
 """
 
 import os
+import unittest
 
 import numpy as np
 
+from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
 from opencsp.app.sofast.lib.ImageCalibrationScaling import ImageCalibrationScaling
 from opencsp.app.sofast.lib.MeasurementSofastFringe import (
     MeasurementSofastFringe as Measurement,
 )
 from opencsp.app.sofast.lib.ProcessSofastFringe import ProcessSofastFringe as Sofast
 from opencsp.common.lib.camera.Camera import Camera
-from opencsp.app.sofast.lib.DisplayShape import DisplayShape as Display
+from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
 from opencsp.common.lib.tool.hdf5_tools import load_hdf5_datasets
 
 
-def test_undefined():
-    # Get test data location
-    base_dir = os.path.join(opencsp_code_dir(), 'test/data/measurements_sofast_fringe')
-
-    # Directory Setup
-    file_dataset = os.path.join(base_dir, 'calculations_undefined_mirror/data.h5')
-    file_measurement = os.path.join(base_dir, 'measurement_facet.h5')
-
-    # Load data
-    camera = Camera.load_from_hdf(file_dataset)
-    display = Display.load_from_hdf(file_dataset)
-    measurement = Measurement.load_from_hdf(file_measurement)
-    calibration = ImageCalibrationScaling.load_from_hdf(file_dataset)
-
-    # Calibrate measurement
-    measurement.calibrate_fringe_images(calibration)
-
-    # Load calculation/user data
-    datasets = [
-        'DataSofastCalculation/facet/facet_000/slopes_facet_xy',
-        'DataSofastCalculation/facet/facet_000/slope_coefs_facet',
-        'DataSofastInput/surface_params/facet_000/initial_focal_lengths_xy',
-        'DataSofastInput/surface_params/facet_000/robust_least_squares',
-        'DataSofastInput/surface_params/facet_000/downsample',
-        'DataSofastInput/surface_params/facet_000/surface_type',
-    ]
-    data = load_hdf5_datasets(datasets, file_dataset)
-
-    # Load sofast params
-    datasets = [
-        'DataSofastInput/sofast_params/mask_hist_thresh',
-        'DataSofastInput/sofast_params/mask_filt_width',
-        'DataSofastInput/sofast_params/mask_filt_thresh',
-        'DataSofastInput/sofast_params/mask_thresh_active_pixels',
-        'DataSofastInput/sofast_params/mask_keep_largest_area',
-        'DataSofastInput/sofast_params/perimeter_refine_axial_search_dist',
-        'DataSofastInput/sofast_params/perimeter_refine_perpendicular_search_dist',
-        'DataSofastInput/sofast_params/facet_corns_refine_step_length',
-        'DataSofastInput/sofast_params/facet_corns_refine_perpendicular_search_dist',
-        'DataSofastInput/sofast_params/facet_corns_refine_frac_keep',
-    ]
-    params = load_hdf5_datasets(datasets, file_dataset)
-
-    # Instantiate sofast object
-    sofast = Sofast(measurement, camera, display)
-
-    # Update parameters
-    sofast.params.mask_hist_thresh = params['mask_hist_thresh']
-    sofast.params.mask_filt_width = params['mask_filt_width']
-    sofast.params.mask_filt_thresh = params['mask_filt_thresh']
-    sofast.params.mask_thresh_active_pixels = params['mask_thresh_active_pixels']
-    sofast.params.mask_keep_largest_area = params['mask_keep_largest_area']
-    sofast.params.geometry_params.perimeter_refine_axial_search_dist = params[
-        'perimeter_refine_axial_search_dist'
-    ]
-    sofast.params.geometry_params.perimeter_refine_perpendicular_search_dist = params[
-        'perimeter_refine_perpendicular_search_dist'
-    ]
-    sofast.params.geometry_params.facet_corns_refine_step_length = params[
-        'facet_corns_refine_step_length'
-    ]
-    sofast.params.geometry_params.facet_corns_refine_perpendicular_search_dist = params[
-        'facet_corns_refine_perpendicular_search_dist'
-    ]
-    sofast.params.geometry_params.facet_corns_refine_frac_keep = params[
-        'facet_corns_refine_frac_keep'
-    ]
-
-    # Define surface data
-    if data['surface_type'] == 'parabolic':
-        surface_data = dict(
-            surface_type=data['surface_type'],
+class test_IntegrationUndefined(unittest.TestCase):
+    def test_undefined(self):
+        # Get test data location
+        base_dir = os.path.join(opencsp_code_dir(), 'test/data/measurements_sofast_fringe')
+
+        # Directory Setup
+        file_dataset = os.path.join(base_dir, 'calculations_undefined_mirror/data.h5')
+        file_measurement = os.path.join(base_dir, 'measurement_facet.h5')
+
+        # Load data
+        camera = Camera.load_from_hdf(file_dataset)
+        display = Display.load_from_hdf(file_dataset)
+        measurement = Measurement.load_from_hdf(file_measurement)
+        calibration = ImageCalibrationScaling.load_from_hdf(file_dataset)
+
+        # Calibrate measurement
+        measurement.calibrate_fringe_images(calibration)
+
+        # Load calculation/user data
+        datasets = [
+            'DataSofastCalculation/facet/facet_000/slopes_facet_xy',
+            'DataSofastCalculation/facet/facet_000/slope_coefs_facet',
+            'DataSofastInput/surface_params/facet_000/initial_focal_lengths_xy',
+            'DataSofastInput/surface_params/facet_000/robust_least_squares',
+            'DataSofastInput/surface_params/facet_000/downsample',
+        ]
+        data = load_hdf5_datasets(datasets, file_dataset)
+
+        # Load sofast params
+        datasets = [
+            'DataSofastInput/sofast_params/mask_hist_thresh',
+            'DataSofastInput/sofast_params/mask_filt_width',
+            'DataSofastInput/sofast_params/mask_filt_thresh',
+            'DataSofastInput/sofast_params/mask_thresh_active_pixels',
+            'DataSofastInput/sofast_params/mask_keep_largest_area',
+            'DataSofastInput/sofast_params/perimeter_refine_axial_search_dist',
+            'DataSofastInput/sofast_params/perimeter_refine_perpendicular_search_dist',
+            'DataSofastInput/sofast_params/facet_corns_refine_step_length',
+            'DataSofastInput/sofast_params/facet_corns_refine_perpendicular_search_dist',
+            'DataSofastInput/sofast_params/facet_corns_refine_frac_keep',
+        ]
+        params = load_hdf5_datasets(datasets, file_dataset)
+
+        # Instantiate sofast object
+        sofast = Sofast(measurement, camera, display)
+
+        # Update parameters
+        sofast.params.mask_hist_thresh = params['mask_hist_thresh']
+        sofast.params.mask_filt_width = params['mask_filt_width']
+        sofast.params.mask_filt_thresh = params['mask_filt_thresh']
+        sofast.params.mask_thresh_active_pixels = params['mask_thresh_active_pixels']
+        sofast.params.mask_keep_largest_area = params['mask_keep_largest_area']
+        sofast.params.geometry_params.perimeter_refine_axial_search_dist = params[
+            'perimeter_refine_axial_search_dist'
+        ]
+        sofast.params.geometry_params.perimeter_refine_perpendicular_search_dist = params[
+            'perimeter_refine_perpendicular_search_dist'
+        ]
+        sofast.params.geometry_params.facet_corns_refine_step_length = params[
+            'facet_corns_refine_step_length'
+        ]
+        sofast.params.geometry_params.facet_corns_refine_perpendicular_search_dist = params[
+            'facet_corns_refine_perpendicular_search_dist'
+        ]
+        sofast.params.geometry_params.facet_corns_refine_frac_keep = params[
+            'facet_corns_refine_frac_keep'
+        ]
+
+        # Define surface data
+        surface = Surface2DParabolic(
             initial_focal_lengths_xy=data['initial_focal_lengths_xy'],
             robust_least_squares=bool(data['robust_least_squares']),
             downsample=data['downsample'],
         )
-    else:
-        surface_data = dict(
-            surface_type=data['surface_type'],
-            robust_least_squares=bool(data['robust_least_squares']),
-            downsample=data['downsample'],
-        )
 
-    # Run SOFAST
-    sofast.process_optic_undefined(surface_data)
+        # Run SOFAST
+        sofast.process_optic_undefined(surface)
 
-    # Test
-    slopes = sofast.data_characterization_facet[0].slopes_facet_xy
-    slope_coefs = sofast.data_characterization_facet[0].slope_coefs_facet
+        # Test
+        slopes = sofast.data_characterization_facet[0].slopes_facet_xy
+        slope_coefs = sofast.data_characterization_facet[0].slope_coefs_facet
 
-    np.testing.assert_allclose(data['slopes_facet_xy'], slopes, atol=1e-7, rtol=0)
-    np.testing.assert_allclose(
-        data['slope_coefs_facet'], slope_coefs, atol=1e-8, rtol=0
-    )
+        np.testing.assert_allclose(data['slopes_facet_xy'], slopes, atol=1e-7, rtol=0)
+        np.testing.assert_allclose(
+            data['slope_coefs_facet'], slope_coefs, atol=1e-8, rtol=0
+        )
 
 
 if __name__ == '__main__':
-    test_undefined()
-    print('All tests run.')
+    unittest.main()
diff --git a/opencsp/app/sofast/test/test_project_fixed_pattern_target.py b/opencsp/app/sofast/test/test_project_fixed_pattern_target.py
index 322ae7c55..112259443 100644
--- a/opencsp/app/sofast/test/test_project_fixed_pattern_target.py
+++ b/opencsp/app/sofast/test/test_project_fixed_pattern_target.py
@@ -3,6 +3,7 @@
 
 import os
 
+import matplotlib
 import pytest
 
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
diff --git a/opencsp/common/lib/deflectometry/ParamsSlopeSolverAbstract.py b/opencsp/common/lib/deflectometry/ParamsSlopeSolverAbstract.py
new file mode 100644
index 000000000..342a5a719
--- /dev/null
+++ b/opencsp/common/lib/deflectometry/ParamsSlopeSolverAbstract.py
@@ -0,0 +1,11 @@
+from abc import ABC
+from dataclasses import dataclass
+
+
+@dataclass
+class ParamsSlopeSolverAbstract(ABC):
+    """Abstract SlopeSolver input parameters class. Contains parameters
+    common to all surface types.
+    """
+    robust_least_squares: bool
+    downsample: int
diff --git a/opencsp/common/lib/deflectometry/ParamsSlopeSolverParaboloid.py b/opencsp/common/lib/deflectometry/ParamsSlopeSolverParaboloid.py
new file mode 100644
index 000000000..2add6c2ef
--- /dev/null
+++ b/opencsp/common/lib/deflectometry/ParamsSlopeSolverParaboloid.py
@@ -0,0 +1,10 @@
+from dataclasses import dataclass
+
+from opencsp.common.lib.deflectometry.ParamsSlopeSolverAbstract import ParamsSlopeSolverAbstract
+
+
+@dataclass
+class ParamsSlopeSolverParaboloid(ParamsSlopeSolverAbstract):
+    """SlopeSolver input parameters class for parabolic surface type
+    """
+    initial_focal_lengths_xy: tuple[float, float]
diff --git a/opencsp/common/lib/deflectometry/ParamsSlopeSolverPlano.py b/opencsp/common/lib/deflectometry/ParamsSlopeSolverPlano.py
new file mode 100644
index 000000000..88e9c9439
--- /dev/null
+++ b/opencsp/common/lib/deflectometry/ParamsSlopeSolverPlano.py
@@ -0,0 +1,9 @@
+from dataclasses import dataclass
+
+from opencsp.common.lib.deflectometry.ParamsSlopeSolverAbstract import ParamsSlopeSolverAbstract
+
+
+@dataclass
+class ParamsSlopeSolverPlano(ParamsSlopeSolverAbstract):
+    """SlopeSolver input parameters class for plano (perfectly flat) surface type
+    """
diff --git a/opencsp/common/lib/deflectometry/SlopeSolver.py b/opencsp/common/lib/deflectometry/SlopeSolver.py
index 2439a9d39..b88579ea6 100644
--- a/opencsp/common/lib/deflectometry/SlopeSolver.py
+++ b/opencsp/common/lib/deflectometry/SlopeSolver.py
@@ -7,14 +7,16 @@
 from opencsp.common.lib.deflectometry.SlopeSolverDataDebug import SlopeSolverDataDebug
 from opencsp.common.lib.deflectometry.SlopeSolverData import SlopeSolverData
 import opencsp.common.lib.deflectometry.slope_fitting_2d as sf2
-from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
-from opencsp.common.lib.deflectometry.Surface2DPlano import Surface2DPlano
+from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract
 from opencsp.common.lib.geometry.Uxyz import Uxyz
 from opencsp.common.lib.geometry.Vxyz import Vxyz
 from opencsp.common.lib.geometry.TransformXYZ import TransformXYZ
 
 
 class SlopeSolver:
+    """Class that solves for the surface slopes of optics in deflectometry
+    systems."""
+
     def __init__(
         self,
         v_optic_cam_optic: Vxyz,
@@ -24,9 +26,9 @@ def __init__(
         v_optic_screen_optic: Vxyz,
         v_align_point_optic: Vxyz,
         dist_optic_screen: float,
-        surface_data: dict,
+        surface: Surface2DAbstract,
         debug: SlopeSolverDataDebug = SlopeSolverDataDebug(),
-    ):
+    ) -> 'SlopeSolver':
         """
         Initializes the slope solving object.
 
@@ -46,62 +48,13 @@ def __init__(
             Position of align point in optic coordinates.
         dist_optic_screen : float
             Measured optic to screen distance.
-        surface_data : dict
-            Dictionary containing surface data information to use when
-            solving slopes. The data fields depend on the surface fit being
-            performed. The following options are supported:
-
-            1) Parabolic fit
-                surface_type
-                initial_focal_lengths_xy
-                robust_least_squares
-                downsample
-            2) Plano fit
-                surface_type
-                robust_least_squares
-                downsample
-            3) Spherical fit (will raise NotImplementedError)
-                surface_type
-                radius
-                robust_least_squares
-                downsample
-
-            Data field descriptions
-            -----------------------
-                surface_type : str {'parabolic', 'plano', 'spherical'}
-                    The type of surface being characterized
-                initial_focal_lengths_xy : list[tuple(float, float)]
-                    The focal lengths to use as the starting point for the
-                    fitting algorithm.
-                robust_least_squares : bool
-                    To use robust least squares fitting, or just least squares
-                    fitting.
-                downsample : int
-                    The amount to downsample data for surface fitting
-                radius : float
-                    The initial radius to use as the starting point for the
-                    fitting algorithm.
+        surface : Surface2DAbstract
+            2D surface definition class.
         debug: SlopeSolverDataDebug
             SlopeSolverDataDebug object for debugging.
-
         """
-        # Instantiate surface fit object depending on fit type
-        surface_data_copy = surface_data.copy()
-        self.surface_type = surface_data_copy.pop('surface_type')
-        if self.surface_type == 'parabolic':
-            self.surface = Surface2DParabolic(**surface_data_copy)
-        elif self.surface_type == 'plano':
-            self.surface = Surface2DPlano(**surface_data_copy)
-        elif self.surface_type == 'spherical':
-            raise NotImplementedError(
-                'Currently, "spherical" surface type is not implemented.'
-            )
-        else:
-            raise ValueError(
-                f'Given surface_type "{self.surface_type:s}" not supported.'
-            )
-
         # Store inputs in class
+        self.surface = surface
         self.v_optic_cam_optic = v_optic_cam_optic
         self.u_active_pixel_pointing_optic = u_active_pixel_pointing_optic
         self.u_measure_pixel_pointing_optic = u_measure_pixel_pointing_optic
diff --git a/opencsp/common/lib/deflectometry/Surface2DAbstract.py b/opencsp/common/lib/deflectometry/Surface2DAbstract.py
index 4218d4e60..34ebc5d90 100644
--- a/opencsp/common/lib/deflectometry/Surface2DAbstract.py
+++ b/opencsp/common/lib/deflectometry/Surface2DAbstract.py
@@ -1,4 +1,4 @@
-from abc import abstractmethod, ABC
+from abc import abstractmethod
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -6,9 +6,10 @@
 
 from opencsp.common.lib.geometry.Uxyz import Uxyz
 from opencsp.common.lib.geometry.Vxyz import Vxyz
+from opencsp.common.lib.tool.hdf5_tools import HDF5_IO_Abstract
 
 
-class Surface2DAbstract(ABC):
+class Surface2DAbstract(HDF5_IO_Abstract):
     """Representation of 2d surface for SOFAST processing"""
 
     def __init__(self):
diff --git a/opencsp/common/lib/deflectometry/Surface2DParabolic.py b/opencsp/common/lib/deflectometry/Surface2DParabolic.py
index 9e44bd9b7..bcb8801b8 100644
--- a/opencsp/common/lib/deflectometry/Surface2DParabolic.py
+++ b/opencsp/common/lib/deflectometry/Surface2DParabolic.py
@@ -5,6 +5,7 @@
 from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract
 from opencsp.common.lib.geometry.Uxyz import Uxyz
 from opencsp.common.lib.geometry.Vxyz import Vxyz
+from opencsp.common.lib.tool.hdf5_tools import save_hdf5_datasets, load_hdf5_datasets
 
 
 class Surface2DParabolic(Surface2DAbstract):
@@ -327,3 +328,35 @@ def shift_all(self, v_align_optic_step: Vxyz) -> None:
         self.v_optic_cam_optic += v_align_optic_step
         self.v_screen_points_optic += v_align_optic_step
         self.v_optic_screen_optic += v_align_optic_step
+
+    def save_to_hdf(self, file: str, prefix: str = ''):
+        data = [
+            self.initial_focal_lengths_xy,
+            self.robust_least_squares,
+            self.downsample,
+            'parabolic'
+        ]
+        datasets = [
+            prefix + 'ParamsSurface/initial_focal_lengths_xy',
+            prefix + 'ParamsSurface/robust_least_squares',
+            prefix + 'ParamsSurface/downsample',
+            prefix + 'ParamsSurface/surface_type',
+        ]
+        save_hdf5_datasets(data, datasets, file)
+
+    @classmethod
+    def load_from_hdf(cls, file: str, prefix: str = ''):
+        # Check surface type
+        data = load_hdf5_datasets([prefix + 'ParamsSurface/surface_type'], file)
+        if data['surface_type'] != 'parabolic':
+            raise ValueError(
+                f'Surface2DParabolic cannot load surface type, {data["surface_type"]:s}')
+
+        # Load
+        datasets = [
+            prefix + 'ParamsSurface/initial_focal_lengths_xy',
+            prefix + 'ParamsSurface/robust_least_squares',
+            prefix + 'ParamsSurface/downsample',
+        ]
+        data = load_hdf5_datasets(datasets, file)
+        return cls(**data)
diff --git a/opencsp/common/lib/deflectometry/Surface2DPlano.py b/opencsp/common/lib/deflectometry/Surface2DPlano.py
index 41d97851d..2a2a00803 100644
--- a/opencsp/common/lib/deflectometry/Surface2DPlano.py
+++ b/opencsp/common/lib/deflectometry/Surface2DPlano.py
@@ -5,6 +5,7 @@
 from opencsp.common.lib.deflectometry.Surface2DAbstract import Surface2DAbstract
 from opencsp.common.lib.geometry.Uxyz import Uxyz
 from opencsp.common.lib.geometry.Vxyz import Vxyz
+from opencsp.common.lib.tool.hdf5_tools import save_hdf5_datasets, load_hdf5_datasets
 
 
 class Surface2DPlano(Surface2DAbstract):
@@ -230,3 +231,32 @@ def shift_all(self, v_align_optic_step: Vxyz) -> None:
         self.u_active_pixel_pointing_optic += v_align_optic_step
         self.u_measure_pixel_pointing_optic += v_align_optic_step
         self.v_optic_screen_optic += v_align_optic_step
+
+    def save_to_hdf(self, file: str, prefix: str = ''):
+        data = [
+            self.robust_least_squares,
+            self.downsample,
+            'parabolic'
+        ]
+        datasets = [
+            prefix + 'ParamsSurface/robust_least_squares',
+            prefix + 'ParamsSurface/downsample',
+            prefix + 'ParamsSurface/surface_type',
+        ]
+        save_hdf5_datasets(data, datasets, file)
+
+    @classmethod
+    def load_from_hdf(cls, file: str, prefix: str = ''):
+        # Check surface type
+        data = load_hdf5_datasets([prefix + 'ParamsSurface/surface_type'], file)
+        if data['surface_type'] != 'parabolic':
+            raise ValueError(
+                f'Surface2DPlano cannot load surface type, {data["surface_type"]:s}')
+
+        # Load
+        datasets = [
+            prefix + 'ParamsSurface/robust_least_squares',
+            prefix + 'ParamsSurface/downsample',
+        ]
+        data = load_hdf5_datasets(datasets, file)
+        return cls(**data)
diff --git a/opencsp/common/lib/deflectometry/test/test_SlopeSolver.py b/opencsp/common/lib/deflectometry/test/test_SlopeSolver.py
index baab9c364..beb5e0ddb 100644
--- a/opencsp/common/lib/deflectometry/test/test_SlopeSolver.py
+++ b/opencsp/common/lib/deflectometry/test/test_SlopeSolver.py
@@ -13,6 +13,7 @@
 )
 from opencsp.app.sofast.lib.SpatialOrientation import SpatialOrientation
 from opencsp.common.lib.deflectometry.SlopeSolver import SlopeSolver
+from opencsp.common.lib.deflectometry.Surface2DParabolic import Surface2DParabolic
 from opencsp.common.lib.geometry.Uxyz import Uxyz
 from opencsp.common.lib.geometry.Vxyz import Vxyz
 from opencsp.common.lib.opencsp_path.opencsp_root_path import opencsp_code_dir
@@ -55,19 +56,11 @@ def setUpClass(cls):
         ori.orient_optic_cam(r_cam_optic, v_cam_optic_cam)
 
         # Perform calculations
-        if data['surface_type'] == 'parabolic':
-            surface_data = dict(
-                surface_type=data['surface_type'],
-                initial_focal_lengths_xy=data['initial_focal_lengths_xy'],
-                robust_least_squares=bool(data['robust_least_squares']),
-                downsample=data['downsample'],
-            )
-        elif data['surface_type'] == 'plano':
-            surface_data = dict(
-                surface_type=data['surface_type'],
-                robust_least_squares=bool(data['robust_least_squares']),
-                downsample=data['downsample'],
-            )
+        surface = Surface2DParabolic(
+            initial_focal_lengths_xy=data['initial_focal_lengths_xy'],
+            robust_least_squares=bool(data['robust_least_squares']),
+            downsample=data['downsample'],
+        )
         kwargs = {
             'v_optic_cam_optic': ori.v_optic_cam_optic,
             'u_active_pixel_pointing_optic': Uxyz(data['u_pixel_pointing_facet']),
@@ -76,7 +69,7 @@ def setUpClass(cls):
             'v_optic_screen_optic': ori.v_optic_screen_optic,
             'v_align_point_optic': measurement.measure_point,
             'dist_optic_screen': measurement.optic_screen_dist,
-            'surface_data': surface_data,
+            'surface': surface,
         }
 
         # Solve slopes
@@ -118,3 +111,7 @@ def test_slopes(self):
         np.testing.assert_allclose(
             data['slopes_facet_xy'], self.data_slope.slopes_facet_xy, atol=1e-8, rtol=0
         )
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/opencsp/common/lib/deflectometry/test/test_Surface2D.py b/opencsp/common/lib/deflectometry/test/test_Surface2D.py
index 1d6bdf6bb..d34c67eb7 100644
--- a/opencsp/common/lib/deflectometry/test/test_Surface2D.py
+++ b/opencsp/common/lib/deflectometry/test/test_Surface2D.py
@@ -1,6 +1,6 @@
 """Unit test suite to test Surface2D type classes
 """
-
+from os.path import dirname, join
 import unittest
 
 import numpy as np
@@ -10,12 +10,17 @@
 from opencsp.common.lib.deflectometry.Surface2DPlano import Surface2DPlano
 from opencsp.common.lib.geometry.Uxyz import Uxyz
 from opencsp.common.lib.geometry.Vxyz import Vxyz
+import opencsp.common.lib.tool.file_tools as ft
 
 
 class Test2DSurface(unittest.TestCase):
     @classmethod
-    def setup_class(cls):
+    def setUpClass(cls):
+        # Generate test data
         cls.data_test = [generate_2DParabolic(), generate_2DPlano()]
+        # Save location
+        cls.dir_save = join(dirname(__file__), 'data/output')
+        ft.create_directories_if_necessary(cls.dir_save)
 
     def test_intersect(self):
         """Tests the intersection of rays with fit surface."""
@@ -93,6 +98,17 @@ def test_fit_normal(self):
                 # Test
                 np.testing.assert_allclose(n_design.data, data_exp.data)
 
+    def test_io(self):
+        """Test saving to HDF5"""
+        prefix = 'test_folder/'
+        for idx, surf in enumerate(self.data_test):
+            surf_cur: Surface2DAbstract = surf[0]
+            file = join(self.dir_save, f'test_surface_{idx:d}.h5')
+            # Test saving
+            surf_cur.save_to_hdf(file, prefix)
+            # Test loading
+            surf_cur.load_from_hdf(file, prefix)
+
 
 def generate_2DParabolic() -> (
     tuple[Surface2DParabolic, Vxyz, np.ndarray, np.ndarray, np.ndarray, Uxyz, Uxyz]
@@ -246,12 +262,4 @@ def generate_2DPlano() -> (
 
 
 if __name__ == '__main__':
-    Test = Test2DSurface()
-    Test.setup_class()
-
-    Test.test_intersect()
-    Test.test_calculate_slopes()
-    Test.test_fit_slopes()
-    Test.test_fit_surf()
-    Test.test_design_normal()
-    Test.test_fit_normal()
+    unittest.main()
diff --git a/opencsp/common/lib/tool/hdf5_tools.py b/opencsp/common/lib/tool/hdf5_tools.py
index 32f9ce0df..389144b80 100644
--- a/opencsp/common/lib/tool/hdf5_tools.py
+++ b/opencsp/common/lib/tool/hdf5_tools.py
@@ -1,6 +1,8 @@
+from abc import abstractmethod, ABC
+import os
+
 import h5py
 import numpy as np
-import os
 
 import opencsp.common.lib.tool.file_tools as ft
 import opencsp.common.lib.tool.image_tools as it
@@ -205,8 +207,7 @@ def unzip(hdf5_path_name_ext: str, destination_dir: str, dataset_format='npy'):
                 aspect_ratio = max(shape[0], shape[1]) / min(shape[0], shape[1])
                 if (shape[0] >= 10 and shape[1] >= 10) and (aspect_ratio < 10.001):
                     dataset_path_name_ext = _create_dataset_path(
-                        hdf5_dir, possible_images[i][0], ".png"
-                    )
+                        hdf5_dir, possible_images[i][0], ".png")
                     # assumed grayscale or RGB
                     if (len(shape) == 2) or (shape[2] in [1, 3]):
                         img = it.numpy_to_image(np_image)
@@ -252,3 +253,36 @@ def unzip(hdf5_path_name_ext: str, destination_dir: str, dataset_format='npy'):
             np.savetxt(dataset_path_name + ".csv", squeezed, delimiter=",")
 
     return hdf5_dir
+
+
+class HDF5_SaveAbstract(ABC):
+    """Abstract class for saving to HDF5 format"""
+
+    @abstractmethod
+    def save_to_hdf(self, file: str, prefix: str = '') -> None:
+        """Saves data to given file. Data is stored as: PREFIX + Folder/Field_1
+
+        Parameters
+        ----------
+        file : str
+            HDF file to save to
+        prefix : str
+            Prefix to append to folder path within HDF file (folders must be separated by "/")
+        """
+
+
+class HDF5_IO_Abstract(HDF5_SaveAbstract):
+    """Abstract class for loading from HDF5 format"""
+
+    @classmethod
+    @abstractmethod
+    def load_from_hdf(cls, file: str, prefix: str = ''):
+        """Loads data from given file. Assumes data is stored as: PREFIX + Folder/Field_1
+
+        Parameters
+        ----------
+        file : str
+            HDF file to save to
+        prefix : str
+            Prefix to append to folder path within HDF file (folders must be separated by "/")
+        """