diff --git a/falco/config/ModelParameters.py b/falco/config/ModelParameters.py
index 3441450..c72d48a 100755
--- a/falco/config/ModelParameters.py
+++ b/falco/config/ModelParameters.py
@@ -9,6 +9,7 @@
 from falco.config.yaml_loader import object_constructor, load_from_str
 from falco.util import _spec_arg
 from falco.config import Probe, ProbeSchedule, Object
+import astropy
 
 
 class ModelParameters(Object):
@@ -251,7 +252,7 @@ def from_yaml(text, context=None):
 
         data = load_from_str(
             text,
-            {'np': numpy, 'falco': falco, 'math': math},
+            {'np': numpy, 'falco': falco, 'math': math, 'astropy': astropy},
             context,
             {
                 "!Probe": object_constructor(Probe),
diff --git a/falco/config/TestbedInterface.py b/falco/config/TestbedInterface.py
new file mode 100644
index 0000000..47ce0b2
--- /dev/null
+++ b/falco/config/TestbedInterface.py
@@ -0,0 +1,9 @@
+from abc import ABC, abstractmethod
+
+class TestbedInterface(ABC):
+    '''A general purpose interface for all testbeds to implement.'''
+
+    @abstractmethod
+    def get_sbp_image(self, subband_index):
+        pass
+
diff --git a/falco/ctrl.py b/falco/ctrl.py
index 5c83f1a..6735630 100644
--- a/falco/ctrl.py
+++ b/falco/ctrl.py
@@ -14,7 +14,7 @@
 import falco
 
 
-def wrapper(mp, cvar, jacStruct):
+def wrapper(mp, cvar, jacStruct, tb=None):
     """
     Outermost wrapper function for all the controller functions.
 
@@ -98,7 +98,7 @@ def wrapper(mp, cvar, jacStruct):
         print('Control beginning ...')
         # Established, conventional controllers
         if mp.controller.lower() == 'plannedefc':
-            dDM = _planned_efc(mp, cvar)
+            dDM = _planned_efc(mp, cvar, tb=tb)
         elif mp.controller.lower() == 'gridsearchefc':
             dDM = _grid_search_efc(mp, cvar)
 
@@ -556,7 +556,7 @@ def _planned_ad_efc(mp, cvar):
     return dDM
 
 
-def _planned_efc(mp, cvar):
+def _planned_efc(mp, cvar, tb=None):
     """
     Perform a scheduled/planned set of EFC iterations.
 
@@ -627,7 +627,7 @@ def _planned_efc(mp, cvar):
 
         for ni in range(Nvals):
 
-            [InormVec[ni], dDM_temp] = _efc(ni, vals_list, mp, cvar)
+            [InormVec[ni], dDM_temp] = _efc(ni, vals_list, mp, cvar, tb=tb)
             ImCube[ni, :, :] = dDM_temp.Itotal
 
             # delta voltage commands
@@ -703,7 +703,7 @@ def _planned_efc(mp, cvar):
         vals_list = [(x, y) for y in np.array([cvar.latestBestDMfac])
                      for x in np.array([log10regSchedOut])]
 
-        [cvar.cMin, dDM] = _efc(ni, vals_list, mp, cvar)
+        [cvar.cMin, dDM] = _efc(ni, vals_list, mp, cvar, tb=tb)
         cvar.Im = np.squeeze(dDM.Itotal)
         if mp.ctrl.flagUseModel:
             print(('Model expects scheduled log10(reg) = %.1f\t to give ' +
@@ -800,7 +800,7 @@ def efc_schedule_generator(scheduleMatrix):
     return Nitr, relinItrVec, gridSearchItrVec, log10regSched, dm_ind_sched
 
 
-def _efc(ni, vals_list, mp, cvar):
+def _efc(ni, vals_list, mp, cvar, tb=None):
     """
     Compute the main EFC equation. Called by a wrapper controller function.
 
@@ -856,7 +856,7 @@ def _efc(ni, vals_list, mp, cvar):
         Itotal = falco.imaging.get_expected_summed_image(mp, cvar, dDM)
         InormAvg = np.mean(Itotal[mp.Fend.corr.maskBool])
     else:  # Get actual image from full model or testbed
-        Itotal = falco.imaging.get_summed_image(mp)
+        Itotal = falco.imaging.get_summed_image(mp, tb=tb)
         InormAvg = np.mean(Itotal[mp.Fend.corr.maskBool])
     dDM.Itotal = Itotal
 
diff --git a/falco/est.py b/falco/est.py
index 6e38f8e..81f7b89 100644
--- a/falco/est.py
+++ b/falco/est.py
@@ -9,7 +9,7 @@
 from . import check
 
 
-def wrapper(mp, ev, jacStruct):
+def wrapper(mp, ev, jacStruct, tb = None):
     """
     Select and run the chosen estimator.
 
@@ -23,6 +23,8 @@ def wrapper(mp, ev, jacStruct):
         Structure containing estimation variables.
     jacStruct : ModelParameters
         Structure containing control Jacobians for each specified DM.
+    tb : falco.config.TestbedInterface or None
+        (Optional) Control interface for a physical testbed.
 
     Returns
     -------
@@ -33,15 +35,15 @@ def wrapper(mp, ev, jacStruct):
 
         falco.est.perfect(mp, ev)
         with falco.util.TicToc('Getting updated summed image'):
-            ev.Im = falco.imaging.get_summed_image(mp)
+            ev.Im = falco.imaging.get_summed_image(mp, tb)
 
     elif mp.estimator in ['pairwise', 'pairwise-square', 'pairwise-rect',
                           'pwp-bp-square', 'pwp-bp', 'pwp-kf']:
 
         if mp.est.flagUseJac:  # Send in the Jacobian if true
-            falco.est.pairwise_probing(mp, ev, jacStruct=jacStruct)
+            falco.est.pairwise_probing(mp, ev, jacStruct=jacStruct, tb=tb)
         else:  # Otherwise don't pass the Jacobian
-            falco.est.pairwise_probing(mp, ev)
+            falco.est.pairwise_probing(mp, ev, tb=tb)
 
     return None
 
@@ -219,7 +221,7 @@ def _est_perfect_Efield_with_Zernikes_in_parallel(mp, ilist, inds_list):
     return E2D
 
 
-def pairwise_probing(mp, ev, jacStruct=np.array([])):
+def pairwise_probing(mp, ev, jacStruct=np.array([]), tb = None):
     """
     Estimate the dark hole E-field with pair-wise probing.
 
@@ -230,6 +232,8 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
     ev : falco.config.Object()
     jacStruct : array_like, optional
         Array containing the control Jacobian. Default is an empty array.
+    tb : falco.config.TestbedInterface or None
+        (Optional) Control interface for a physical testbed.
 
     Returns
     -------
@@ -394,7 +398,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
 
             # Take initial, unprobed image (for unprobed DM settings).
             whichImage = 0
-            I0 = falco.imaging.get_sbp_image(mp, iSubband)
+            I0 = falco.imaging.get_sbp_image(mp, iSubband, tb)
             I0vec = I0[mp.Fend.corr.maskBool]  # Vectorize the dark hole
 
             # Image already includes all stars, so don't sum over star loop
@@ -460,7 +464,7 @@ def pairwise_probing(mp, ev, jacStruct=np.array([])):
                     mp.dm2.V = DM2Vnom + dDM2Vprobe
 
                 # Take probed image
-                Im = falco.imaging.get_sbp_image(mp, iSubband)
+                Im = falco.imaging.get_sbp_image(mp, iSubband, tb)
 
                 # ImNonneg = Im
                 # ImNonneg[Im < 0] = 0
@@ -875,4 +879,4 @@ def gen_pairwise_probe(mp, InormDes, phaseShift, starIndex, rotation):
     probeCmd = falco.dm.fit_surf_to_act(dm, probeHeight)
     probeCmd = mp.est.probe.gainFudge[starIndex] * probeCmd  # Scale the probe amplitude empirically if needed
 
-    return probeCmd
\ No newline at end of file
+    return probeCmd
diff --git a/falco/imaging.py b/falco/imaging.py
index 316435c..d8b49f8 100644
--- a/falco/imaging.py
+++ b/falco/imaging.py
@@ -246,7 +246,7 @@ def _model_full_norm_wrapper(mp, ilist, inds_list):
     return np.max(np.abs(Etemp)**2)
 
 
-def get_summed_image(mp):
+def get_summed_image(mp, tb = None):
     """
     Get the broadband image over the entire bandpass.
 
@@ -257,6 +257,8 @@ def get_summed_image(mp):
     ----------
     mp: falco.config.ModelParameters
         Structure of model parameters
+    tb: falco.config.TestbedInterface
+        (Optional) control interface for a physical testbed
 
     Returns
     -------
@@ -270,7 +272,7 @@ def get_summed_image(mp):
     if not (mp.flagParallel and mp.flagSim):
         summedImage = 0
         for si in range(mp.Nsbp):
-            summedImage += mp.sbp_weights[si] * get_sbp_image(mp, si)
+            summedImage += mp.sbp_weights[si] * get_sbp_image(mp, si, tb)
 
     else:  # Compute simulated images in parallel
 
@@ -360,7 +362,7 @@ def _get_single_sim_full_image(mp, ilist, vals_list):
     return np.abs(Estar)**2  # Apply spectral weighting outside this function
 
 
-def get_sbp_image(mp, si):
+def get_sbp_image(mp, si, tb = None):
     """
     Get an image in the specified sub-bandpass.
 
@@ -373,6 +375,8 @@ def get_sbp_image(mp, si):
         Structure of model parameters
     si: int
         Index of sub-bandpass for which to take the image
+    tb: falco.config.TestbedInterface or None
+        (Optional) Control interface for a physical testbed
 
     Returns
     -------
@@ -387,7 +391,8 @@ def get_sbp_image(mp, si):
     if mp.flagSim:
         Isbp = get_sim_sbp_image(mp, si)
     else:
-        Isbp = get_testbed_sbp_image(mp, si)
+        tb.dm.apply(mp.dm1.V)
+        Isbp = tb.get_sbp_image(si)
 
     return Isbp
 
diff --git a/falco/wfsc.py b/falco/wfsc.py
index cae6fc2..625f0c5 100644
--- a/falco/wfsc.py
+++ b/falco/wfsc.py
@@ -9,7 +9,7 @@
 import falco
 
 
-def loop(mp, out):
+def loop(mp, out, tb = None):
     """
     Loop over the estimator and controller for WFSC.
 
@@ -19,6 +19,8 @@ def loop(mp, out):
         Structure of model parameters
     out : falco.config.Object
         Output variables
+    tb : falco.config.TestbedInterface or None
+        (Optional) Control interface for a physical testbed.
 
     Returns
     -------
@@ -90,7 +92,7 @@ def loop(mp, out):
         if Itr > 0:
             EestPrev = ev.Eest  # save previous estimate for Delta E plot
 
-        falco.est.wrapper(mp, ev, jacStruct)
+        falco.est.wrapper(mp, ev, jacStruct, tb)
 
         store_intensities(mp, out, ev, Itr)
 
@@ -122,7 +124,7 @@ def loop(mp, out):
 
         # Send a copy of jacStruct so that spatial weights don't show up 
         # outside the controller or get applied multiple times.
-        falco.ctrl.wrapper(mp, cvar, copy(jacStruct))
+        falco.ctrl.wrapper(mp, cvar, copy(jacStruct), tb=tb)
 
         # Store key data in out object
         out.log10regHist[Itr] = cvar.log10regUsed