diff --git a/optvl/optvl_class.py b/optvl/optvl_class.py
index 24d55eb2..4a7cc3e0 100644
--- a/optvl/optvl_class.py
+++ b/optvl/optvl_class.py
@@ -216,6 +216,11 @@ class OVLSolver(object):
         "XYZref":  ["CASE_R","XYZREF0"],
         "CDp":   ["CASE_R", "CDREF0"],
     }
+    
+    viz_aifoil_data_list = [
+        "casec",
+        "tasec"
+    ]
 
     # fmt: on
 
@@ -230,7 +235,7 @@ def __init__(
         debug: Optional[bool] = False,
         timing: Optional[bool] = False,
     ):
-        """Initalize the python and fortran libary from the given objects
+        """Initialize the python and fortran library from the given objects
 
         Args:
             geo_file: AVL geometry file
@@ -326,10 +331,13 @@ def __init__(
 
         # control surfaces added in __init__
         # TODO: the keys of this dict aren't used
-        self.con_var_to_fort_var = {
-            "alpha": ["CASE_R", "ALFA"],
-            "beta": ["CASE_R", "BETA"],
-        }
+        self.con_var_list = [
+            "alpha",
+            "beta",
+            "roll rate",
+            "pitch rate",
+            "yaw rate",
+        ]
 
         control_names = self.get_control_names()
         self.dindex_to_con_surf = OrderedDict()
@@ -342,7 +350,7 @@ def __init__(
         idx_control_start = np.max([x for x in self.conval_idx_dict.values()]) + 1
         for idx_c_var, c_name in enumerate(control_names):
             self.conval_idx_dict[c_name] = idx_control_start + idx_c_var
-            self.con_var_to_fort_var[c_name] = ["CASE_R", "DELCON"]
+            self.con_var_list.append(c_name)
 
         var_to_suffix = {
             "alpha": "AL",
@@ -467,8 +475,8 @@ def _init_map_data(self):
 
     def _setup_surface_maps(self, surf_name: str, idx_surf: int, num_sec: int):
         """Used by the init_map_data and load_input_dict functions to generate which slices of the Fortran array for a
-        given geometry, panneling, control surface, or design variable correspond to the given surface. This data is
-        stored the surf_geom_to_fort_var dictionary.
+        given geometry, paneling, control surface, or design variable correspond to the given surface. This data is
+        stored in the surf_geom_to_fort_var dictionary.
 
         Args:
             surf_name: The name of the surface
@@ -562,7 +570,7 @@ def _setup_surface_maps(self, surf_name: str, idx_surf: int, num_sec: int):
 
     def _setup_body_maps(self, body_name: str, idx_body: int):
         """Used by the init_map_data and load_input_dict functions to generate which slices of the Fortran array for a
-        given geometry or discretization variable correspond to the given body. This data is stored the
+        given geometry or discretization variable correspond to the given body. This data is stored in the
         body_param_to_fort_var dictionary.
 
         Args:
@@ -582,7 +590,7 @@ def _setup_body_maps(self, body_name: str, idx_body: int):
 
     def _setup_section_maps(self, surf_name: str, idx_surf: int, num_sec: int, nasec_arr: np.ndarray):
         """Used by the init_map_data and load_input_dict functions to generate which slices of the Fortran array for a
-        given section geometry variable correspond to the given surface and section. This data is stored the
+        given section geometry variable correspond to the given surface and section. This data is stored in the
         surf_section_geom_to_fort_var dictionary.
 
         Args:
@@ -628,7 +636,7 @@ def load_input_dict(self, input_dict: dict, pre_check: bool = True, post_check:
         Args:
             input_dict: input dictionary in optvl format
             pre_check: perform additional verification of the user's input dictionary before loading into AVL
-            post_check: verify certain inputs values are correctly reflected in the Fortran layer
+            post_check: verify certain input values are correctly reflected in the Fortran layer
         """
 
         # Initialize Variables and Counters
@@ -942,6 +950,7 @@ def check_type(key, avl_vars, given_val, cast_type=True):
 
                         for key in self.surf_section_geom_to_fort_var[surf_name]:
                             avl_vars_secs = self.surf_section_geom_to_fort_var[surf_name][key]
+                            
                             avl_vars = (avl_vars_secs[0], avl_vars_secs[1], avl_vars_secs[2][j])
 
                             if key not in surf_dict:
@@ -1172,7 +1181,7 @@ def set_mesh(self, idx_surf: int, mesh: np.ndarray, flatten:bool=True, update_nv
         Args:
             idx_surf (int): the surface to apply the mesh to
             mesh (np.ndarray): XYZ mesh array (nx,ny,3)
-            flatten (bool): Should OptVL flatten the mesh when placing vorticies and control points
+            flatten (bool): Should OptVL flatten the mesh when placing vortices and control points
             update_nvs (bool): Should OptVL update the number of spanwise elements for the given mesh
             update_nvc (bool): Should OptVL update the number of chordwise elements for the given mesh
         """
@@ -1208,7 +1217,6 @@ def set_mesh(self, idx_surf: int, mesh: np.ndarray, flatten:bool=True, update_nv
 
     def get_mesh(self, idx_surf: int, concat_dup_mesh: bool = False):
         """Returns the current set mesh coordinates from AVL as a numpy array.
-        Note this is intended for 
 
         Args:
             idx_surf (int): the surface to get the mesh for
@@ -1260,7 +1268,7 @@ def set_section_naca(self, isec: int, isurf: int, nasec: int, naca: str, xfminma
             isec: section number to set the airfoil mesh
             isurf: surface number to set the airfoil mesh
             nasec: number of points to evaluate the interpolated camber line and thickness curves at
-            naca: 4-digit naca specificaion as a string
+            naca: 4-digit naca specification as a string
             xfminmax: length 2 array with the min and max x/c to slice the airfoil
         """
 
@@ -1311,7 +1319,7 @@ def set_section_coordinates(
             isurf: surface number to set the airfoil mesh
             nasec: number of points to evaluate the interpolated camber line and thickness curves at
             x: airfoil x-coordinate array
-            y: airfoil y-coodinate array
+            y: airfoil y-coordinate array
             xfminmax: length 2 array with the min and max x/c to slice the airfoil
             storecoords: store the raw input coordinates in common block
         """
@@ -1335,10 +1343,10 @@ def set_body_coordinates(self, ibod: int, nasec: int, x: np.ndarray, y: np.ndarr
 
 
         Args:
-            ibod: body number to set the outer mold line too
+            ibod: body number to set the outer mold line to
             nasec: number of points to evaluate the interpolated camber line and thickness curves at
             x: oml x-coordinate array
-            y: oml y-coodinate array
+            y: oml y-coordinate array
             xfminmax: length 2 array with the min and max x/c to slice the oml
             storecoords: store the raw input coordinates in common block
         """
@@ -1478,16 +1486,17 @@ def execute_run(self, tol: float = 0.00002):
         """Run the analysis (equivalent to the AVL command `x` in the OPER menu)
 
         Args:
-            tol: the tolerace of the Newton solver used for triming the aircraft
+            tol: the tolerance of the Newton solver used for trimming the aircraft
         """
         self.set_avl_fort_arr("CASE_R", "EXEC_TOL", tol)
         self.avl.oper()
 
     def set_variable(self, var: str, val: float):
-        """set a variable for the run case (equivalent to setting a variable in AVL's OPER menu)
+        """Set a variable for the run case (equivalent to setting a variable in AVL's OPER menu)
+
         Args:
-            var: variable to be constrained ["alpha"", "beta"", "roll rate", "pitch rate", "yaw rate"] or any control surface.
-            val: target value of `con_var`
+            var: variable to be set ["alpha", "beta", "roll rate", "pitch rate", "yaw rate"] or any control surface.
+            val: target value of `var`
         """
         avl_variables = {
             "alpha": ("CASE_R", "ALFA"),
@@ -1532,11 +1541,12 @@ def set_variable(self, var: str, val: float):
         # set the contraint value so that the set value is used in analysis
         self.set_avl_fort_arr("CASE_R", "CONVAL", val, (self.conval_idx_dict[var],))
 
-    def get_variable(self, var: str):
-        """set a variable for the run case (equivalent to setting a variable in AVL's OPER menu)
+    def get_variable(self, var: str, inRadians: bool = False):
+        """Get a variable for the run case (equivalent to reading a variable in AVL's OPER menu)
+
         Args:
-            var: variable to be constrained ["alpha"", "beta"", "roll rate", "pitch rate", "yaw rate"] or any control surface.
-            val: target value of `con_var`
+            var: variable to retrieve ["alpha", "beta", "roll rate", "pitch rate", "yaw rate"] or any control surface.
+            inRadians: if True, return alpha and beta in radians instead of degrees. Roll, pitch, and yaw rates are always dimensionless.
         """
         avl_variables = {
             "alpha": ("CASE_R", "ALFA"),
@@ -1559,6 +1569,10 @@ def get_variable(self, var: str):
             "yaw rate": 2 / bref,
         }
 
+        if inRadians:
+            vars_factors["alpha"] = 1.0
+            vars_factors["beta"] = 1.0
+
         # Check that the variable is valid
         if var in avl_variables:
             # save the name of the avl_var
@@ -1579,9 +1593,9 @@ def set_constraint(self, var: str, con_var: str, val: float):
         """Set the constraints on the analysis case (equivalent to setting a constraint in AVL's OPER menu)
 
         Args:
-            var: variable to be constrained ["alpha"", "beta"", "roll rate", "pitch rate", "yaw rate"] or any control surface.
+            var: variable to be constrained ["alpha", "beta", "roll rate", "pitch rate", "yaw rate"] or any control surface.
             val: target value of `con_var`
-            con_var: variable output that needs to be constrained. It could be any value for `var` plus ["CL", "CY", "Cl", "Cm", "Cn"]. If None, than `var` is also the `con_var`
+            con_var: variable output that needs to be constrained. It could be any value for `var` plus ["CL", "CY", "Cl", "Cm", "Cn"]. If None, then `var` is also the `con_var`
 
         """
         avl_variables = {
@@ -1635,7 +1649,7 @@ def set_constraint(self, var: str, con_var: str, val: float):
         self.avl.conset(avl_var, f"{avl_con_var} {val} \n")
 
     def set_trim_condition(self, variable: str, val: float):
-        """Set a variable of the trim condition (analogus to the AVL's C1 command from the OPER menu)
+        """Set a variable of the trim condition (analogous to the AVL's C1 command from the OPER menu)
 
         Args:
             variable: variable to be set. Options are ["bankAng", "CL", "velocity", "mass", "dens", "G", "X cg","Y cg","Z cg"]
@@ -1666,7 +1680,7 @@ def get_total_forces(self) -> Dict[str, float]:
         """Get the aerodynamic data for the last run case and return it as a dictionary.
 
         Returns:
-            Dict[str, float]: Dictionary of aerodynamic data. The keys the aerodyanmic coefficients.
+            Dict[str, float]: Dictionary of aerodynamic data. The keys are the aerodynamic coefficients.
         """
 
         total_data = {}
@@ -1681,10 +1695,10 @@ def get_total_forces(self) -> Dict[str, float]:
         return total_data
 
     def get_body_forces(self) -> Dict[str, float]:
-        """Get the aerodynamic data for the last run case and return it as a dictionary.
+        """Get the aerodynamic force data for each body from the last run case.
 
         Returns:
-            Dict[str, float]: Dictionary of aerodynamic data. The keys the aerodyanmic coefficients.
+            Dict[str, float]: Dictionary of aerodynamic data. The keys are the aerodynamic coefficients.
         """
 
         body_data = {}
@@ -1708,7 +1722,7 @@ def get_control_stab_derivs(self) -> Dict[str, float]:
         for the current analysis run
 
         Returns:
-            stab_deriv_dict: The dictionary of control surface derivatives, d{force coefficent}/d{control surface}.
+            stab_deriv_dict: The dictionary of control surface derivatives, d{force coefficient}/d{control surface}.
         """
 
         deriv_data = {}
@@ -1725,7 +1739,7 @@ def get_control_stab_derivs(self) -> Dict[str, float]:
         return deriv_data
 
     def get_stab_derivs(self) -> Dict[str, float]:
-        """gets the stability derivates after an analysis run
+        """Gets the stability derivatives after an analysis run
 
         Returns:
             stab_deriv_dict: Dictionary of stability derivatives.
@@ -1739,10 +1753,10 @@ def get_stab_derivs(self) -> Dict[str, float]:
         return deriv_data
 
     def get_body_axis_derivs(self) -> Dict[str, float]:
-        """gets the body-axis derivates after an analysis run
+        """Gets the body-axis derivatives after an analysis run
 
         Returns:
-            body_deriv_dict: Dictionary of stability derivatives.
+            body_deriv_dict: Dictionary of body-axis derivatives.
         """
         deriv_data = {}
 
@@ -1773,11 +1787,11 @@ def set_reference_data(self, ref_data: Dict[str, float], set_default_value:bool
         return ref_data
 
     def get_avl_fort_arr(self, common_block: str, variable: str, slicer: Optional[slice] = None) -> np.ndarray:
-        """Get data from the Fortran level common block data structure. see AVL.INC for all availible variables
+        """Get data from the Fortran level common block data structure. See AVL.INC for all available variables
 
         Args:
             common_block: Name of the common block of the variable like `CASE_R`
-            variable: Name of the variable to retrive
+            variable: Name of the variable to retrieve
             slicer: slice applied to the common block variable to return a subset of the data. i.e. (100) or slice(2, 5)
 
         Returns:
@@ -1806,11 +1820,11 @@ def get_avl_fort_arr(self, common_block: str, variable: str, slicer: Optional[sl
         return val
 
     def set_avl_fort_arr(self, common_block: str, variable: str, val: float, slicer: Optional[slice] = None) -> None:
-        """Set data from the Fortran level common block data structure. see AVL.INC for all availible variables
+        """Set data from the Fortran level common block data structure. See AVL.INC for all available variables
 
         Args:
             common_block: Name of the common block of the variable like `CASE_R`
-            variable: Name of the variable to retrive
+            variable: Name of the variable to retrieve
             val: value to set, which can be a numpy array
             slicer: slice applied to the common block variable to return a subset of the data. i.e. (100) or slice(2, 5)
 
@@ -1853,7 +1867,7 @@ def set_avl_fort_arr(self, common_block: str, variable: str, val: float, slicer:
         return
 
     def get_surface_forces(self) -> Dict[str, Dict[str, float]]:
-        """Returns the force data from each surface (including mirriored surfaces)
+        """Returns the force data from each surface (including mirrored surfaces)
 
         Returns:
             surf_data_dict: a dictionary of surface data where the first key is the surface and the second is the force coefficient
@@ -1881,7 +1895,7 @@ def get_surface_forces(self) -> Dict[str, Dict[str, float]]:
 
     def get_parameter(self, param_key: str) -> float:
         """
-        Analogous to ruinont Modify parameters for the OPER menu to view parameters.
+        Analogous to the Modify parameters (M) command from the OPER menu in AVL, but for reading parameters.
 
         Args:
             param_key: the name of the parameter to return
@@ -1970,7 +1984,7 @@ def get_control_deflections(self) -> Dict[str, float]:
         return def_dict
 
     def get_control_deflection(self, con_surf) -> Dict[str, float]:
-        """get the deflections of the control surfaces
+        """Get the deflection of the specified control surface
 
         Returns:
             val: deflection of control surface
@@ -1985,10 +1999,11 @@ def get_control_deflection(self, con_surf) -> Dict[str, float]:
         return val
 
     def set_control_deflection(self, con_surf, val) -> Dict[str, float]:
-        """set the deflections of the control surfaces
+        """Set the deflection of a control surface
 
-        args:
-            con_surf : Name or D value (D1, D2, etc.) of the control surface
+        Args:
+            con_surf: Name or D value (D1, D2, etc.) of the control surface
+            val: deflection value to set
         """
         con_surf_names = list(self.con_surf_to_dindex.keys())
         con_surf_dindex = list(self.dindex_to_con_surf.keys())
@@ -2005,9 +2020,9 @@ def set_control_deflection(self, con_surf, val) -> Dict[str, float]:
         self.set_avl_fort_arr("CASE_R", "DELCON", val, slicer=(idx_con_surf,))
 
     def set_control_deflections(self, def_dict: Dict[str, float]):
-        """get the deflections of all the control surfaces
+        """Set the deflections of all the control surfaces
 
-        args:
+        Args:
             def_dict: dictionary of control surfaces as the keys and deflections as the values
         """
         control_surfaces = self.get_control_names()
@@ -2020,7 +2035,7 @@ def get_hinge_moments(self) -> Dict[str, float]:
         """Get the hinge moments from the fortran layer and return them as a dictionary
 
         Returns:
-            hinge_moments: array of control surface moments. The order the control surfaces are declared are the indices,
+            hinge_moments: dictionary of control surface hinge moments keyed by control surface name
         """
         hinge_moments = {}
 
@@ -2173,7 +2188,7 @@ def get_eigenvalues(self) -> np.ndarray:
         return eig_vals
 
     def get_eigenvectors(self) -> np.ndarray:
-        """After running an eigenmode calculation, this function will return the eigenvalues in the order used by AVL
+        """After running an eigenmode calculation, this function will return the eigenvectors in the order used by AVL
 
         Returns:
             eig_vec: 2D array of eigen vectors
@@ -2199,13 +2214,13 @@ def get_system_matrix(self, in_body_axis=False) -> np.ndarray:
         return Asys
 
     def get_system_matrices(self, in_body_axis=False) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
-        """returns the A,B, and R system matrices used in amode.f
-        
-        args:
+        """Returns the A, B, and R system matrices used in amode.f
+
+        Args:
             in_body_axis: apply the sign changes to the matrix to put it in the body axis
-        
+
         Returns:
-            Asys: 2D array representing the system matrix for the eigen value analysis
+            Asys: 2D array representing the system matrix for the eigenvalue analysis
             Bsys: 2D array representing the system matrix for control surfaces
             Rsys: 1D array representing the RHS of the dynamics equation
         """
@@ -2284,10 +2299,10 @@ def get_control_names(self) -> List[str]:
         return control_names
 
     def get_design_var_names(self) -> List[str]:
-        """Get the names of the design_var surfaces
+        """Get the names of the design variables
 
         Returns:
-            design_var_names: list of design_var surface names
+            design_var_names: list of design variable names
         """
         fort_names = self.get_avl_fort_arr("CASE_C", "GNAME")
         design_var_names = self._fort_char_array_to_str_list(fort_names)
@@ -2297,7 +2312,7 @@ def get_surface_names(self, remove_dublicated: Optional[bool] = False) -> List[s
         """Get the surface names from the geometry
 
         Args:
-            remove_dublicated: remove the surface that were created by duplication about symmetry planes
+            remove_dublicated: remove the surfaces that were created by duplication about symmetry planes
 
         Returns:
             surf_names: list of surface names
@@ -2321,7 +2336,7 @@ def get_body_names(self, remove_dublicated: Optional[bool] = False) -> List[str]
         """Get the body names from the geometry
 
         Args:
-            remove_dublicated: remove the body that were created by duplication about symmetry planes
+            remove_dublicated: remove the bodies that were created by duplication about symmetry planes
 
         Returns:
             body_names: list of body names
@@ -2368,12 +2383,12 @@ def get_con_surf_param(self, surf_name: str, idx_sec: int, param: str) -> np.nda
         return param
 
     def __get_des_var_param(self, surf_name: str, idx_sec: int, param: str) -> np.ndarray:
-        """Returns the parameters that define the control surface. Can also get design variables (AVL).
+        """Returns the parameters that define a design variable.
 
         Args:
-            surf_name: the name of the surface containing the control surface
-            idx_sec: the section index of the control surface data
-            param: control surface parameter to get
+            surf_name: the name of the surface containing the design variable
+            idx_sec: the section index of the design variable data
+            param: design variable parameter to get
 
         Returns:
             parm: parameter value
@@ -2464,12 +2479,12 @@ def get_surface_param(self, surf_name: str, param: str) -> np.ndarray:
         return copy.deepcopy(param)  # return the value of the array, but not a reference to avoid sideffects
 
     def set_surface_param(self, surf_name: str, param: str, val: float, update_geom: bool = True):
-        """Set a parameter of a specified surface. Supports setting params related to geometry and panelling.
+        """Set a parameter of a specified surface. Supports setting params related to geometry and paneling.
         Section geometry can be directly set here but this is not recommended. Use set_section_coordinates instead.
 
         Args:
             surf_name: the surface containing the parameter
-            param: the surface parameter to return. Could be either geometric or paneling
+            param: the surface parameter to set. Could be either geometric or paneling
             val: value to set
             update_geom: flag to update the geometry after setting
         """
@@ -2539,7 +2554,7 @@ def get_surface_params(
         include_des_vars: bool = False,
         include_airfoils: bool = False,
     ) -> Dict[str, Dict[str, Any]]:
-        """Get all the surface level parameters for each suface
+        """Get all the surface level parameters for each surface
 
         Args:
             include_geom: flag to include geometry data in the output. The data is ["scale", "translate", "angle", "xles", "yles", "zles", "chords", "aincs", "clcdsec", "claf"]
@@ -2548,7 +2563,7 @@ def get_surface_params(
             include_con_surf:  flag to include control surface and design variable data in the output. This is data like the hinge vector and gain.
             include_airfoils:  flag to include airfoil file data in the output
 
-        Return:
+        Returns:
             surf_data: Nested dictionary where the 1st key is the surface name and the 2nd key is the parameter.
         """
         surf_data = {}
@@ -2686,7 +2701,7 @@ def get_body_param(self, body_name: str, param: str) -> np.ndarray:
 
 
         Returns:
-            val: the val of parameter of the body
+            val: the value of the parameter of the body
         """
 
         body_names = self.get_body_names()
@@ -2720,7 +2735,7 @@ def set_body_param(self, body_name: str, param: str, val, update_geom: bool = Tr
 
         Args:
             body_name: the body containing the parameter
-            param: the surface parameter to return. Could be either geometric or paneling
+            param: the body parameter to set. Could be either geometric or paneling
             val: value to set
             update_geom: flag to update the geometry after setting
         """
@@ -2802,7 +2817,7 @@ def get_body_params(
         return body_data
 
     def set_body_params(self, body_data: Dict[str, Dict[str, Any]]):
-        """Set the give body data of the current geometry.
+        """Set the given body data of the current geometry.
 
         Args:
             body_data: Nested dictionary where the 1st key is the body name and the 2nd key is the parameter.
@@ -3135,12 +3150,11 @@ def _write_data(key_list: List[str], newline: bool = True):
 
     # region --- Utility functions
     def get_num_surfaces(self) -> int:
-        """Returns the number of surface including duplicated
+        """Returns the number of surfaces including duplicated surfaces
 
         Returns:
             val: number of surfaces
         """
-        """Get the number of surfaces in the geometry"""
         return self.get_avl_fort_arr("CASE_I", "NSURF")
 
     def get_surface_index(self, surf_name: str) -> int:
@@ -3176,7 +3190,7 @@ def get_num_sections(self, surf_name: str) -> int:
             surf_name: name of the surface
 
         Returns:
-            nsec: numer of sections
+            nsec: number of sections
         """
         idx_surf = self.get_surface_index(surf_name)
         slice_idx_surf = (idx_surf,)
@@ -3206,7 +3220,7 @@ def get_mesh_size(self) -> int:
         return int(self.get_avl_fort_arr("CASE_I", "NVOR"))
 
     def get_mesh_data(self) -> Dict[str, int]:
-        """Get the number of vortices in the mesh
+        """Get mesh size data (bodies, surfaces, strips, vortices, and variables)
 
         Returns:
             mesh_data : dict
@@ -3224,8 +3238,8 @@ def get_mesh_data(self) -> Dict[str, int]:
         return mesh_data
 
     def _str_to_fort_str(self, py_string, num_max_char):
-        """Setting arrays of strings in Fortran can be kinda nasty. This
-        takes a strings and returns the char array.
+        """Setting arrays of strings in Fortran is complex. This
+        takes a string and returns the char array.
         """
 
         arr = np.zeros((), dtype=f"|S{num_max_char}")
@@ -3239,8 +3253,8 @@ def _str_to_fort_str(self, py_string, num_max_char):
         return arr
 
     def _str_to_fort_char_array(self, py_string, num_max_char):
-        """Setting arrays of strings in Fortran can be kinda nasty. This
-        takes a strings and returns the char array.
+        """Setting arrays of strings in Fortran is complex. This
+        takes a string and returns the char array.
         """
 
         arr = np.zeros(1, dtype=f"|S{num_max_char}")
@@ -3254,8 +3268,8 @@ def _str_to_fort_char_array(self, py_string, num_max_char):
         return arr
 
     def _str_list_to_fort_char_array(self, strList, num_max_char):
-        """Setting arrays of strings in Fortran can be kinda nasty. This
-        takes a list of strings and returns the array.
+        """Setting arrays of strings in Fortran is complex. This
+        takes a list of strings and returns the char array.
         """
 
         arr = np.zeros((len(strList), num_max_char), dtype="str")
@@ -3326,11 +3340,39 @@ def _get_deriv_key(self, var: str, func: str) -> str:
     # ---------------------------
     # --- Derivative routines ---
     # ---------------------------
+    
+    # --- utils ---
+    def zero_seed_dict(self, seed_dict: Dict) -> Dict:
+        for key, value in seed_dict.items():
+            seed_dict[key] = self._zero_value(value)
+        return seed_dict
+
+    def _zero_value(self, value):
+        if isinstance(value, dict):
+            for k, v in value.items():
+                value[k] = self._zero_value(v)
+            return value
+        elif isinstance(value, list):
+            return [self._zero_value(v) for v in value]
+        elif isinstance(value, np.ndarray):
+            value[...] = 0
+            return value
+        else:
+            return 0
+    
+    def deep_update(self, d, other):
+        for k, v in other.items():
+            if isinstance(v, dict) and isinstance(d.get(k), dict):
+                self.deep_update(d[k], v)
+            else:
+                d[k] = v
+        return d
+    
 
     # --- input ad seeds ---
     def get_variable_ad_seeds(self) -> Dict[str, float]:
         var_seeds = {}
-        for con in self.con_var_to_fort_var:
+        for con in self.con_var_list:
             idx_con = self.conval_idx_dict[con]
             blk = "CASE_R" + self.ad_suffix
             var = "CONVAL" + self.ad_suffix
@@ -3360,7 +3402,7 @@ def set_variable_ad_seeds(self, con_seeds: Dict[str, Dict[str, float]], mode: st
                 self.set_avl_fort_arr(blk, var, val, slicer=slicer)
 
             elif mode == "FD":
-                # # reverse lookup in the con_var_to_fort_var dict
+                # # reverse lookup in the con_var_list dict
                 if con in self.con_surf_to_dindex:
                     val = self.get_control_deflection(con)
                 else:
@@ -3454,23 +3496,56 @@ def get_geom_ad_seeds(self) -> Dict[str, Dict[str, float]]:
                 var += self.ad_suffix
 
                 geom_seeds[surf_key][geom_key] = copy.deepcopy(self.get_avl_fort_arr(blk, var, slicer=slicer))
+                
+            for geom_key in self.surf_section_geom_to_fort_var[surf_key]:
+                if geom_key in self.viz_aifoil_data_list:
+                    continue
+                
+                blk, var, sec_slicers = self.surf_section_geom_to_fort_var[surf_key][geom_key]
+                blk += self.ad_suffix
+                var += self.ad_suffix
+                
+                geom_seeds[surf_key][geom_key] = []
+                for sec_slice in (sec_slicers):
+                    geom_seeds[surf_key][geom_key].append(copy.deepcopy(self.get_avl_fort_arr(blk, var, slicer=sec_slice)))
 
         return geom_seeds
 
     def set_geom_ad_seeds(self, geom_seeds: Dict[str, float], mode: str = "AD", scale=1.0) -> None:
         for surf_key in geom_seeds:
             for geom_key in geom_seeds[surf_key]:
-                blk, var, slicer = self.surf_geom_to_fort_var[surf_key][geom_key]
+                if geom_key in self.surf_geom_to_fort_var[surf_key]:
+                    blk, var, slicer = self.surf_geom_to_fort_var[surf_key][geom_key]
+                    
+                    if mode == "AD":
+                        blk += self.ad_suffix
+                        var += self.ad_suffix
+                        val = geom_seeds[surf_key][geom_key] * scale
+                    elif mode == "FD":
+                        val = self.get_avl_fort_arr(blk, var, slicer=slicer)
+                        val += geom_seeds[surf_key][geom_key] * scale
+                    
+                    self.set_avl_fort_arr(blk, var, val, slicer=slicer)
+                
+                elif geom_key in self.viz_aifoil_data_list:
+                    raise ValueError(f"Can not set key {geom_key}. it is only used for viz")
+                elif geom_key in self.surf_section_geom_to_fort_var[surf_key]:
+                    
+                    blk, var, sec_slicers = self.surf_section_geom_to_fort_var[surf_key][geom_key]
 
-                if mode == "AD":
-                    blk += self.ad_suffix
-                    var += self.ad_suffix
-                    val = geom_seeds[surf_key][geom_key] * scale
-                elif mode == "FD":
-                    val = self.get_avl_fort_arr(blk, var, slicer=slicer)
-                    val += geom_seeds[surf_key][geom_key] * scale
-                # print(blk, var, val, slicer)
-                self.set_avl_fort_arr(blk, var, val, slicer=slicer)
+                    if mode == "AD":
+                        blk += self.ad_suffix
+                        var += self.ad_suffix
+                
+                    for idx_slice, sec_slice in enumerate(sec_slicers):
+                        if mode == "AD":
+                            val = geom_seeds[surf_key][geom_key][idx_slice] * scale
+                        elif mode == "FD":
+                            val = self.get_avl_fort_arr(blk, var, slicer=sec_slice)
+                            val += geom_seeds[surf_key][geom_key][idx_slice] * scale
+                
+                        self.set_avl_fort_arr(blk, var, val, slicer=sec_slice)
+                    
 
     def get_mesh_ad_seeds(self) -> Dict[str, Dict[str, float]]:
         mesh_seeds = {}
@@ -3834,7 +3909,7 @@ def _execute_jac_vec_prod_fwd(
             step: Step size to use for the FD mode
 
         Returns:
-            func_seeds: force coifficent AD seeds
+            func_seeds: force coefficient AD seeds
             res_seeds:  residual AD seeds
             consurf_derivs_seeds: Control surface derivatives AD seeds
             stab_derivs_seeds: Stability derivatives AD seeds
@@ -3848,16 +3923,41 @@ def _execute_jac_vec_prod_fwd(
 
         mesh_size = self.get_mesh_size()
         num_control_surfs = self.get_num_control_surfs()
-
-        if con_seeds is None:
-            con_seeds = {}
-
-        if geom_seeds is None:
-            geom_seeds = {}
-
-        if mesh_seeds is None:
-            mesh_seeds = {}
-
+        
+        # start from zero'd con and geom seeds 
+        
+        # --- con seeds ---
+        con_seeds_full = self.get_variable_ad_seeds()
+        self.zero_seed_dict(con_seeds_full)
+        if con_seeds is not None:
+            self.deep_update(con_seeds_full, con_seeds)
+        
+        # --- geom ---
+        geom_seeds_full = self.get_geom_ad_seeds()
+        self.zero_seed_dict(geom_seeds_full)
+        if geom_seeds is not None:
+            self.deep_update(geom_seeds_full, geom_seeds)
+        
+        
+        # --- mesh ---
+        mesh_seeds_full = self.get_mesh_ad_seeds()
+        self.zero_seed_dict(mesh_seeds_full)
+        if mesh_seeds is not None:
+            self.deep_update(mesh_seeds_full, mesh_seeds)
+        
+        # --- param ---
+        param_seeds_full = self.get_parameter_ad_seeds()
+        self.zero_seed_dict(param_seeds_full)
+        if param_seeds is not None:
+            self.deep_update(param_seeds_full, param_seeds)
+        
+        # --- ref ---
+        ref_seeds_full = self.get_reference_ad_seeds()
+        self.zero_seed_dict(ref_seeds_full)
+        if ref_seeds is not None:
+            self.deep_update(ref_seeds_full, ref_seeds)
+        
+        # The arrays can be set to zero if they don't exsist
         if gamma_seeds is None:
             gamma_seeds = np.zeros(mesh_size)
 
@@ -3867,28 +3967,27 @@ def _execute_jac_vec_prod_fwd(
         if gamma_u_seeds is None:
             gamma_u_seeds = np.zeros((self.NUMAX, mesh_size))
 
-        if param_seeds is None:
-            param_seeds = {}
-
-        if ref_seeds is None:
-            ref_seeds = {}
 
         res_slice = (slice(0, mesh_size),)
         res_d_slice = (slice(0, num_control_surfs), slice(0, mesh_size))
         res_u_slice = (slice(0, self.NUMAX), slice(0, mesh_size))
-
+        
         if mode == "AD":
+            # The following body seeds are in the process of being supported
+            # because they are not explicitly passed, zero them out for now
+            body_seeds = np.zeros((self.NLMAX, 3))
+            self.set_avl_fort_arr("VRTX_R_DIFF", "RL_DIFF", body_seeds)
+            
             # set derivative seeds
-            # self.clear_ad_seeds()
-            self.set_variable_ad_seeds(con_seeds)
-            self.set_geom_ad_seeds(geom_seeds)
-            self.set_mesh_ad_seeds(mesh_seeds)
+            self.set_variable_ad_seeds(con_seeds_full)
+            self.set_geom_ad_seeds(geom_seeds_full)
+            self.set_mesh_ad_seeds(mesh_seeds_full)
             self.set_gamma_ad_seeds(gamma_seeds)
             self.set_gamma_d_ad_seeds(gamma_d_seeds)
             self.set_gamma_u_ad_seeds(gamma_u_seeds)
-            self.set_parameter_ad_seeds(param_seeds)
-            self.set_reference_ad_seeds(ref_seeds)
-
+            self.set_parameter_ad_seeds(param_seeds_full)
+            self.set_reference_ad_seeds(ref_seeds_full)
+            
             self.avl.update_surfaces_d()
             self.avl.get_res_d()
             self.avl.velsum_d()
@@ -3903,27 +4002,29 @@ def _execute_jac_vec_prod_fwd(
             res_d_seeds = self.get_residual_d_ad_seeds()
             res_u_seeds = self.get_residual_u_ad_seeds()
 
-            self.set_variable_ad_seeds(con_seeds, scale=0.0)
-            self.set_geom_ad_seeds(geom_seeds, scale=0.0)
-            self.set_mesh_ad_seeds(mesh_seeds, scale=0.0)
+            self.set_variable_ad_seeds(con_seeds_full, scale=0.0)
+            self.set_geom_ad_seeds(geom_seeds_full, scale=0.0)
+            self.set_mesh_ad_seeds(mesh_seeds_full, scale=0.0)
             self.set_gamma_ad_seeds(gamma_seeds, scale=0.0)
             self.set_gamma_d_ad_seeds(gamma_d_seeds, scale=0.0)
             self.set_gamma_u_ad_seeds(gamma_u_seeds, scale=0.0)
-            self.set_parameter_ad_seeds(param_seeds, scale=0.0)
-            self.set_reference_ad_seeds(ref_seeds, scale=0.0)
+            self.set_parameter_ad_seeds(param_seeds_full, scale=0.0)
+            self.set_reference_ad_seeds(ref_seeds_full, scale=0.0)
 
-            # TODO: remove??
+            # TODO: revome this line and see if any of the tests break
+            #  one should not have to zero the gam seed since they are set directly, right?
+            # also the above line should do the zeroing right, set_gamma_ad_seeds
             self.set_avl_fort_arr("VRTX_R_DIFF", "GAM_DIFF", gamma_seeds * 0.0, slicer=res_slice)
 
         if mode == "FD":
-            self.set_variable_ad_seeds(con_seeds, mode="FD", scale=step)
-            self.set_geom_ad_seeds(geom_seeds, mode="FD", scale=step)
-            self.set_mesh_ad_seeds(mesh_seeds, mode="FD", scale=step)
+            self.set_variable_ad_seeds(con_seeds_full, mode="FD", scale=step)
+            self.set_geom_ad_seeds(geom_seeds_full, mode="FD", scale=step)
+            self.set_mesh_ad_seeds(mesh_seeds_full, mode="FD", scale=step)
             self.set_gamma_ad_seeds(gamma_seeds, mode="FD", scale=step)
             self.set_gamma_d_ad_seeds(gamma_d_seeds, mode="FD", scale=step)
             self.set_gamma_u_ad_seeds(gamma_u_seeds, mode="FD", scale=step)
-            self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step)
-            self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step)
+            self.set_parameter_ad_seeds(param_seeds_full, mode="FD", scale=step)
+            self.set_reference_ad_seeds(ref_seeds_full, mode="FD", scale=step)
 
             # propogate the seeds through without resolving
             self.avl.update_surfaces()
@@ -3940,14 +4041,14 @@ def _execute_jac_vec_prod_fwd(
             res_d_peturbed = copy.deepcopy(self.get_avl_fort_arr("VRTX_R", "RES_D", slicer=res_d_slice))
             res_u_peturbed = copy.deepcopy(self.get_avl_fort_arr("VRTX_R", "RES_U", slicer=res_u_slice))
 
-            self.set_variable_ad_seeds(con_seeds, mode="FD", scale=-1 * step)
-            self.set_geom_ad_seeds(geom_seeds, mode="FD", scale=-1 * step)
-            self.set_mesh_ad_seeds(mesh_seeds, mode="FD", scale=-1 * step)
+            self.set_variable_ad_seeds(con_seeds_full, mode="FD", scale=-1 * step)
+            self.set_geom_ad_seeds(geom_seeds_full, mode="FD", scale=-1 * step)
+            self.set_mesh_ad_seeds(mesh_seeds_full, mode="FD", scale=-1 * step)
             self.set_gamma_ad_seeds(gamma_seeds, mode="FD", scale=-1 * step)
             self.set_gamma_d_ad_seeds(gamma_d_seeds, mode="FD", scale=-1 * step)
             self.set_gamma_u_ad_seeds(gamma_u_seeds, mode="FD", scale=-1 * step)
-            self.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step)
-            self.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step)
+            self.set_parameter_ad_seeds(param_seeds_full, mode="FD", scale=-1 * step)
+            self.set_reference_ad_seeds(ref_seeds_full, mode="FD", scale=-1 * step)
 
             self.avl.update_surfaces()
             self.avl.get_res()
@@ -4129,6 +4230,188 @@ def _execute_jac_vec_prod_rev(
 
         return con_seeds, geom_seeds, mesh_seeds, gamma_seeds, gamma_d_seeds, gamma_u_seeds, param_seeds, ref_seeds
 
+    def execute_run_sensitivities_direct(
+        self,
+        con_dvs: Optional[List[str]] = None,
+        geom_dvs: Optional[List[Tuple[(str, str)]]] = None,
+        param_dvs: Optional[List[str]] = None,
+        ref_dvs: Optional[List[str]] = None,
+        add_stab_derivs: Optional[bool] = False,
+        add_body_axis_derivs: Optional[bool] = False,
+        add_consurf_derivs: Optional[bool] = False,
+        print_timings: Optional[bool] = False,
+    ) -> Dict[str, Dict[str, float]]:
+        """Run the sensitivities of the input functionals in adjoint mode
+
+        Returns:
+            sens: a nested dictionary of sensitivities. The first key is the function and the next keys are for the design variables.
+        """
+        
+        sens = {}
+
+        if self.get_avl_fort_arr("CASE_L", "LTIMING"):
+            print_timings = True
+        
+
+        def make_unit_dicts(surfaces, surf_name, key):
+            """Yield dicts with each scalar element set to 1.0, others zero."""
+            val = surfaces[surf_name][key]
+
+            if np.isscalar(val):
+                yield {surf_name: {key: 1.0}}
+
+            elif isinstance(val, np.ndarray):
+                for idx in np.ndindex(val.shape):
+                    out = np.zeros_like(val, dtype=float)
+                    out[idx] = 1.0
+                    yield {surf_name: {key: out}}
+
+            elif isinstance(val, list):
+                # list of arrays (like xasec, sasec)
+                for i, arr in enumerate(val):
+                    for idx in np.ndindex(arr.shape):
+                        out = [np.zeros_like(a, dtype=float) for a in val]
+                        out[i][idx] = 1.0
+                        yield {surf_name: {key: out}}
+
+        # convenience function to get the correct seeds:
+        def get_dv_seed_list():
+            seeds = {}
+            
+            con_seed_list = []
+            if con_dvs is not None:
+                for con in con_dvs:
+                    con_seed_list.append({con:1.0})
+            seeds['con_seeds'] = con_seed_list
+  
+            geom_seed_list = []
+            if geom_dvs is not None:
+                full_geom_seeds = self.get_geom_ad_seeds()
+                for dv in geom_dvs:
+                    for dv_dict in make_unit_dicts(full_geom_seeds, dv[0], dv[1]):
+                        geom_seed_list.append(dv_dict)
+            seeds['geom_seeds'] = geom_seed_list
+            
+            param_seed_list = []
+            if param_dvs is not None:
+                for param in param_dvs:
+                    param_seed_list.append({param:1.0})
+            seeds['param_seeds'] = param_seed_list
+
+            ref_seed_list = []
+            if ref_dvs is not None:
+                for ref in ref_dvs:
+                    ref_seed_list.append({ref:1.0})
+            seeds['ref_seeds'] = ref_seed_list
+                    
+            # for dv in geom_dvs:
+            #     full_seed = full_geom_seeds[dv[0]][dv[1]]
+            #     if isinstance(full_seed, np.ndarray):
+                    
+            #         for i in range(full_seed.size):
+                    
+                
+            return seeds
+            
+            
+        
+        seeds = get_dv_seed_list()
+            
+        
+        def add_deriv_to_dict(dv_seed, sens_dict, seeds):
+            
+            def add_sub_dicts(in_dict, out_dict):
+                # recurse until we don't reach a dictionary
+                for key in in_dict:
+                    val = in_dict[key]
+            
+                    if isinstance(val, dict):
+                        if key not in out_dict:
+                            out_dict[key] = {}
+                        add_sub_dicts(val, out_dict[key])
+                    elif isinstance(val, np.ndarray):
+                      if key not in out_dict:
+                          out_dict[key] = np.zeros_like(val,dtype=float)
+                      out_dict[key] += val * seeds[func_key]
+                    elif isinstance(val, list):
+                        if key not in out_dict:
+                            out_dict[key] = [np.zeros_like(a,dtype=float) for a in val]
+                        for i, arr in enumerate(val):
+                            out_dict[key][i] += arr *seeds[func_key]
+                    else:  # scalar
+                        out_dict[key] = out_dict.get(key, 0.0) + val * seeds[func_key]
+            
+            for func_key in seeds:
+                
+                if func_key not in sens_dict:
+                    sens_dict[func_key] = {}
+                    
+                add_sub_dicts(dv_seed, sens_dict[func_key])
+
+        for dv_type in seeds:
+            for dv_seed in seeds[dv_type]:
+                print(f'processing {dv_seed}')
+                
+                time_last = time.time()
+                # compute dR/dX
+                _, pRpX, _, _, _, pR_dpX, pR_upX = self._execute_jac_vec_prod_fwd(
+                    **{dv_type: dv_seed}
+                )
+                if print_timings:
+                    print(f"Time to get RHS: {time.time() - time_last}")
+                    time_last = time.time()
+
+                # now solve the direct equation
+                # the RHS is absorbing the negative of the total deriv equation
+                # DONT forget the negative
+                self.set_residual_ad_seeds(-1 * pRpX)
+                
+                if add_stab_derivs or add_body_axis_derivs:
+                    self.set_residual_u_ad_seeds(-1 * pR_upX)
+                    solve_res_u_drt = True
+                else:
+                    solve_res_u_drt = False
+                    
+                if add_consurf_derivs:
+                    self.set_residual_d_ad_seeds(-1 * pR_dpX)
+                    solve_res_d_drt = True
+                else:
+                    solve_res_d_drt = False
+
+                # solve only the state direct method
+                self.avl.solve_direct(solve_res_u_drt, solve_res_d_drt)
+                
+                if print_timings:
+                    print(f"Time to solve direct: {time.time() - time_last}")
+                    time_last = time.time()
+
+                # get the resulting adjoint vector (dU/dX) from fortran
+                dUdX = self.get_gamma_ad_seeds()
+                
+                # this is harmless even if we aren't solving for extra derivatives
+                dU_ddX = self.get_gamma_d_ad_seeds()
+                dU_udX = self.get_gamma_u_ad_seeds()
+                
+                # do the last partial derivative to comput the totals with the direct vector
+                func_seeds, _, consurf_derivs_seeds, stab_derivs_seeds, body_axis_derivs_seeds, _, _, = self._execute_jac_vec_prod_fwd(
+                        gamma_seeds=dUdX, gamma_d_seeds=dU_ddX, gamma_u_seeds=dU_udX, **{dv_type: dv_seed}
+                    )
+            
+                
+                add_deriv_to_dict(dv_seed, sens, func_seeds)
+                
+                if add_consurf_derivs:
+                    add_deriv_to_dict(dv_seed, sens, consurf_derivs_seeds)
+                
+                if add_stab_derivs:
+                    add_deriv_to_dict(dv_seed, sens, stab_derivs_seeds)
+                
+                if add_body_axis_derivs:
+                    add_deriv_to_dict(dv_seed, sens, body_axis_derivs_seeds)
+            
+        return sens
+        
+        
     def execute_run_sensitivities(
         self,
         funcs: List[str],
@@ -4143,7 +4426,7 @@ def execute_run_sensitivities(
             funcs: force coefficients to compute the sensitivities with respect to
             stab_derivs: stability derivatives to compute the sensitivities with respect to
             body_axis_derivs: body axis derivatives to compute the sensitivities with respect to
-            consurf_derivs: control surface derivates to compute the sensitivities with respect to
+            consurf_derivs: control surface derivatives to compute the sensitivities with respect to
             print_timings: flag to print timing information
 
         Returns:
@@ -4157,16 +4440,14 @@ def execute_run_sensitivities(
         # set up and solve the adjoint for each function
         for func in funcs:
             sens[func] = {}
-            # get the RHS of the adjoint equation (pFpU)
-            # TODO: remove seeds if it doesn't effect accuracy
-            # self.clear_ad_seeds()
             time_last = time.time()
+            
+            # get the RHS of the adjoint equation (pFpU)
             _, _, _, pfpU, _, _, _, _ = self._execute_jac_vec_prod_rev(func_seeds={func: 1.0})
             if print_timings:
                 print(f"Time to get RHS: {time.time() - time_last}")
                 time_last = time.time()
 
-            # self.clear_ad_seeds()
             # u solver adjoint equation with RHS
             self.set_gamma_ad_seeds(-1 * pfpU)
             solve_gamma_u_adj = False
@@ -4570,7 +4851,7 @@ def plot_geom(self, axes=None, body_color="magenta"):
         """Generate a matplotlib plot of geometry
 
         Args:
-            axes: Matplotlib axis object to add the plots too. If none are given, the axes will be generated.
+            axes: Matplotlib axis object to add the plots to. If none are given, the axes will be generated.
             body_color: Color to use for plotting body geometry (default: 'magenta')
         """
 
@@ -4781,12 +5062,12 @@ def plot_geom_3d(self, axes=None, plot_avl_mesh = True, plot_direct_mesh = False
         By default the flat version of the mesh that satisfies AVL's VLM assumptions is plotted.
         This is either the mesh that comes as a result of AVL's standard geometry specification system
         or the custom user assigned mesh after it has undergone the transformation needed to flatten it to
-        satify the VLM assumptions. There is also an option to overlay the directly assigned mesh. This will
+        satisfy the VLM assumptions. There is also an option to overlay the directly assigned mesh. This will
         plot user assigned mesh as is with no modifications.
 
         Args:
-            axes: Matplotlib axis object to add the plots too. If none are given, the axes will be generated.
-            plot_avl_mesh: If True the AVL flattenned mesh is plotted on the axis
+            axes: Matplotlib axis object to add the plots to. If none are given, the axes will be generated.
+            plot_avl_mesh: If True the AVL flattened mesh is plotted on the axis
             plot_direct_mesh: If True the user assigned mesh will be plotted as is
         """
 
diff --git a/src/aoper.f b/src/aoper.f
index 53f238f9..38685ca5 100644
--- a/src/aoper.f
+++ b/src/aoper.f
@@ -2353,4 +2353,46 @@ subroutine solve_adjoint(solve_stab_deriv_adj, solve_con_surf_adj)
       enddo
       endif
       
+      end !subroutine solve_adjoint
+
+      subroutine solve_direct(solve_stab_deriv_drt, solve_con_surf_drt)
+      ! solves the direct equation with the jacobian to compute dx/du
+      use avl_heap_inc
+      use avl_heap_diff_inc
+      include "AVL.INC"
+      include "AVL_ad_seeds.inc"
+      integer i 
+      logical :: solve_stab_deriv_drt, solve_con_surf_drt
+      
+      CALL SETUP
+      IF(.NOT.LAIC) THEN
+            call factor_AIC
+      ENDIF
+      
+      ! assume we have run the other routines to put dr/dx is in gam_diff
+      do i =1,NVOR
+            GAM_diff(i) = RES_diff(i)
+      enddo
+
+      CALL BAKSUB(NVOR,NVOR,AICN_LU,IAPIV,GAM_diff)
+      
+      if (solve_con_surf_drt) then 
+      DO IC = 1, NCONTROL
+            do i =1,NVOR
+                  GAM_D_diff(i,IC) = RES_D_diff(i,IC)
+            enddo
+            CALL BAKSUB(NVOR,NVOR,AICN_LU,IAPIV,GAM_D_diff(:,IC))
+      enddo
+      endif 
+
+      if (solve_stab_deriv_drt) then 
+      DO IU = 1,6
+            do i =1,NVOR
+                  GAM_U_diff(i,IU) = RES_U_diff(i,IU)
+            enddo
+            
+            CALL BAKSUB(NVOR,NVOR,AICN_LU,IAPIV,GAM_U_diff(:,IU))
+      enddo
+      endif
+      
       end !subroutine solve_adjoint
\ No newline at end of file
diff --git a/src/f2py/libavl.pyf b/src/f2py/libavl.pyf
index ce5919af..67244940 100644
--- a/src/f2py/libavl.pyf
+++ b/src/f2py/libavl.pyf
@@ -117,6 +117,10 @@ python module libavl ! in
         subroutine solve_adjoint(solve_stab_deriv_adj, solve_con_surf_adj)
             logical :: solve_stab_deriv_adj, solve_con_surf_adj
         end subroutine solve_adjoint
+
+        subroutine solve_direct(solve_stab_deriv_drt, solve_con_surf_drt)
+            logical :: solve_stab_deriv_drt, solve_con_surf_drt
+        end subroutine solve_direct
         
         subroutine cpoml(save_file)
             logical :: save_file
diff --git a/tests/test_total_derivs.py b/tests/test_total_derivs.py
index 5d562ada..19892a76 100644
--- a/tests/test_total_derivs.py
+++ b/tests/test_total_derivs.py
@@ -28,11 +28,11 @@ class TestTotals(unittest.TestCase):
     # TODO: beta derivatives likely wrong
 
     def setUp(self):
-        self.ovl_solver = OVLSolver(geo_file=geom_file)
-        self.ovl_solver.set_variable("alpha", 5.0)
-        self.ovl_solver.set_variable("beta", 0.0)
-        self.ovl_solver.set_parameter("Mach", 0.8)
-        self.ovl_solver.execute_run()
+        self.ovl = OVLSolver(geo_file=geom_file)
+        self.ovl.set_variable("alpha", 5.0)
+        self.ovl.set_variable("beta", 0.0)
+        self.ovl.set_parameter("Mach", 0.8)
+        self.ovl.execute_run()
 
     def tearDown(self):
         # Get the memory usage of the current process using psutil
@@ -45,42 +45,42 @@ def finite_dif(self, con_list, geom_seeds, param_seeds, ref_seeds, step=1e-7):
 
         for con in con_list:
             con_seeds[con] = 1.0
-        self.ovl_solver.set_variable_ad_seeds(con_seeds, mode="FD", scale=step)
-        self.ovl_solver.set_geom_ad_seeds(geom_seeds, mode="FD", scale=step)
-        self.ovl_solver.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step)
-        self.ovl_solver.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step)
-
-        self.ovl_solver.avl.update_surfaces()
-        self.ovl_solver.avl.get_res()
-        self.ovl_solver.avl.exec_rhs()
-        self.ovl_solver.avl.get_res()
-        self.ovl_solver.avl.velsum()
-        self.ovl_solver.avl.aero()
-        # self.ovl_solver.execute_run()
-        coef_data_peturb = self.ovl_solver.get_total_forces()
-        consurf_derivs_peturb = self.ovl_solver.get_control_stab_derivs()
-        stab_deriv_derivs_peturb = self.ovl_solver.get_stab_derivs()
-        body_axis_deriv_petrub = self.ovl_solver.get_body_axis_derivs()
-        body_forces_peturb = self.ovl_solver.get_body_forces()
-
-        self.ovl_solver.set_variable_ad_seeds(con_seeds, mode="FD", scale=-1 * step)
-        self.ovl_solver.set_geom_ad_seeds(geom_seeds, mode="FD", scale=-1 * step)
-        self.ovl_solver.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step)
-        self.ovl_solver.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step)
-
-        self.ovl_solver.avl.update_surfaces()
-        self.ovl_solver.avl.get_res()
-        self.ovl_solver.avl.exec_rhs()
-        self.ovl_solver.avl.get_res()
-        self.ovl_solver.avl.velsum()
-        self.ovl_solver.avl.aero()
-        # self.ovl_solver.execute_run()
-
-        coef_data = self.ovl_solver.get_total_forces()
-        consurf_derivs = self.ovl_solver.get_control_stab_derivs()
-        stab_deriv_derivs = self.ovl_solver.get_stab_derivs()
-        body_axis_deriv = self.ovl_solver.get_body_axis_derivs()
-        body_forces = self.ovl_solver.get_body_forces()
+        self.ovl.set_variable_ad_seeds(con_seeds, mode="FD", scale=step)
+        self.ovl.set_geom_ad_seeds(geom_seeds, mode="FD", scale=step)
+        self.ovl.set_parameter_ad_seeds(param_seeds, mode="FD", scale=step)
+        self.ovl.set_reference_ad_seeds(ref_seeds, mode="FD", scale=step)
+
+        self.ovl.avl.update_surfaces()
+        self.ovl.avl.get_res()
+        self.ovl.avl.exec_rhs()
+        self.ovl.avl.get_res()
+        self.ovl.avl.velsum()
+        self.ovl.avl.aero()
+        # self.ovl.execute_run()
+        coef_data_peturb = self.ovl.get_total_forces()
+        consurf_derivs_peturb = self.ovl.get_control_stab_derivs()
+        stab_deriv_derivs_peturb = self.ovl.get_stab_derivs()
+        body_axis_deriv_petrub = self.ovl.get_body_axis_derivs()
+        body_forces_peturb = self.ovl.get_body_forces()
+
+        self.ovl.set_variable_ad_seeds(con_seeds, mode="FD", scale=-1 * step)
+        self.ovl.set_geom_ad_seeds(geom_seeds, mode="FD", scale=-1 * step)
+        self.ovl.set_parameter_ad_seeds(param_seeds, mode="FD", scale=-1 * step)
+        self.ovl.set_reference_ad_seeds(ref_seeds, mode="FD", scale=-1 * step)
+
+        self.ovl.avl.update_surfaces()
+        self.ovl.avl.get_res()
+        self.ovl.avl.exec_rhs()
+        self.ovl.avl.get_res()
+        self.ovl.avl.velsum()
+        self.ovl.avl.aero()
+        # self.ovl.execute_run()
+
+        coef_data = self.ovl.get_total_forces()
+        consurf_derivs = self.ovl.get_control_stab_derivs()
+        stab_deriv_derivs = self.ovl.get_stab_derivs()
+        body_axis_deriv = self.ovl.get_body_axis_derivs()
+        body_forces = self.ovl.get_body_forces()
         
         body_func_seeds = {}
         for body in body_forces:
@@ -111,14 +111,14 @@ def finite_dif(self, con_list, geom_seeds, param_seeds, ref_seeds, step=1e-7):
 
     def test_aero_constraint(self):
         # compare the analytical gradients with finite difference for each constraint and function
-        func_vars = self.ovl_solver.case_var_to_fort_var
-        stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var
-        body_axis_derivs = self.ovl_solver.case_body_derivs_to_fort_var
-        sens_funcs = self.ovl_solver.execute_run_sensitivities(func_vars)
-        sens_sd = self.ovl_solver.execute_run_sensitivities([], stab_derivs=stab_derivs, print_timings=False)
-        sens_bd = self.ovl_solver.execute_run_sensitivities([], body_axis_derivs=body_axis_derivs, print_timings=False)
-
-        for con_key in self.ovl_solver.con_var_to_fort_var:
+        func_vars = self.ovl.case_var_to_fort_var
+        stab_derivs = self.ovl.case_stab_derivs_to_fort_var
+        body_axis_derivs = self.ovl.case_body_derivs_to_fort_var
+        sens_funcs = self.ovl.execute_run_sensitivities(func_vars)
+        sens_sd = self.ovl.execute_run_sensitivities([], stab_derivs=stab_derivs, print_timings=False)
+        sens_bd = self.ovl.execute_run_sensitivities([], body_axis_derivs=body_axis_derivs, print_timings=False)
+
+        for con_key in self.ovl.con_var_to_fort_var:
             # for con_key in ['beta']:
             func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif(
                 [con_key], {}, {}, {}, step=1.0e-5
@@ -209,19 +209,19 @@ def test_geom(self):
         # compare the analytical gradients with finite difference for each
         # geometric variable and function
 
-        surf_key = list(self.ovl_solver.surf_geom_to_fort_var.keys())[0]
-        geom_vars = self.ovl_solver.surf_geom_to_fort_var[surf_key]
-        cs_names = self.ovl_solver.get_control_names()
+        surf_key = list(self.ovl.surf_geom_to_fort_var.keys())[0]
+        geom_vars = self.ovl.surf_geom_to_fort_var[surf_key]
+        cs_names = self.ovl.get_control_names()
 
         consurf_vars = []
-        for func_key in self.ovl_solver.case_derivs_to_fort_var:
-            consurf_vars.append(self.ovl_solver._get_deriv_key(cs_names[0], func_key))
+        for func_key in self.ovl.case_derivs_to_fort_var:
+            consurf_vars.append(self.ovl._get_deriv_key(cs_names[0], func_key))
 
-        func_vars = self.ovl_solver.case_var_to_fort_var
-        stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var
-        body_axis_derivs = self.ovl_solver.case_body_derivs_to_fort_var
+        func_vars = self.ovl.case_var_to_fort_var
+        stab_derivs = self.ovl.case_stab_derivs_to_fort_var
+        body_axis_derivs = self.ovl.case_body_derivs_to_fort_var
 
-        sens = self.ovl_solver.execute_run_sensitivities(
+        sens = self.ovl.execute_run_sensitivities(
             func_vars,
             consurf_derivs=consurf_vars,
             stab_derivs=stab_derivs,
@@ -229,12 +229,12 @@ def test_geom(self):
             print_timings=False,
         )
 
-        # for con_key in self.ovl_solver.con_var_to_fort_var:
+        # for con_key in self.ovl.con_var_to_fort_var:
         sens_FD = {}
-        for surf_key in self.ovl_solver.surf_geom_to_fort_var:
+        for surf_key in self.ovl.surf_geom_to_fort_var:
             sens_FD[surf_key] = {}
             for geom_key in geom_vars:
-                arr = self.ovl_solver.get_surface_param(surf_key, geom_key)
+                arr = self.ovl.get_surface_param(surf_key, geom_key)
                 np.random.seed(arr.size)
                 rand_arr = np.random.rand(*arr.shape)
                 rand_arr /= np.linalg.norm(rand_arr)
@@ -352,12 +352,12 @@ def test_geom(self):
 
     def test_params(self):
         # compare the analytical gradients with finite difference for each constraint and function
-        func_vars = self.ovl_solver.case_var_to_fort_var
-        stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var
+        func_vars = self.ovl.case_var_to_fort_var
+        stab_derivs = self.ovl.case_stab_derivs_to_fort_var
 
-        sens = self.ovl_solver.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs)
+        sens = self.ovl.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs)
 
-        for param_key in self.ovl_solver.param_idx_dict:
+        for param_key in self.ovl.param_idx_dict:
             func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif(
                 [], {}, {param_key: 1.0}, {}, step=1.0e-6
             )
@@ -414,12 +414,12 @@ def test_params(self):
 
     def test_ref(self):
         # compare the analytical gradients with finite difference for each constraint and function
-        func_vars = self.ovl_solver.case_var_to_fort_var
-        stab_derivs = self.ovl_solver.case_stab_derivs_to_fort_var
+        func_vars = self.ovl.case_var_to_fort_var
+        stab_derivs = self.ovl.case_stab_derivs_to_fort_var
 
-        sens = self.ovl_solver.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs)
+        sens = self.ovl.execute_run_sensitivities(func_vars, stab_derivs=stab_derivs)
 
-        for ref_key in self.ovl_solver.ref_var_to_fort_var:
+        for ref_key in self.ovl.ref_var_to_fort_var:
             # for con_key in ['beta']:
             func_seeds, consurf_deriv_seeds, stab_derivs_seeds, body_axis_derivs_seeds = self.finite_dif(
                 [], {}, {}, {ref_key: 1.0}, step=1.0e-5
@@ -477,6 +477,51 @@ def test_ref(self):
                         err_msg=f"{func_key}  wrt {ref_key}",
                     )
 
+class TestDirectVsAdjoint(unittest.TestCase):
+
+    def setUp(self):
+        self.ovl = OVLSolver(geo_file=geom_file)
+        self.ovl.set_variable("alpha", 5.0)
+        self.ovl.set_variable("beta", 0.0)
+        self.ovl.set_parameter("Mach", 0.8)
+        self.ovl.execute_run()
+
+    def tearDown(self):
+        # Get the memory usage of the current process using psutil
+        process = psutil.Process()
+        mb_memory = process.memory_info().rss / (1024 * 1024)  # Convert bytes to MB
+        print(f"{self.id()} Memory usage: {mb_memory:.2f} MB")
+
+    def test_aero_constraint(self):
+        # compare the analytical gradients with finite difference for each constraint and function
+        # func_vars = self.ovl.case_var_to_fort_var
+        stab_derivs = ["dCL/dalpha"]
+        # body_axis_derivs = self.ovl.case_body_derivs_to_fort_var
+        funcs = ["CL", "CD", "Cm"]
+        # funcs = ["CL"]
+        con_dvs = ["alpha"]
+        ref_dvs = ["Sref"]
+        param_dvs = ["Mach"]
+        # geom_dvs = [("Wing", "xles"), ("Wing", "chords")]
+        geom_dvs = [("Wing", "scale"),("Wing", "chords")]
+        sens_adjoint = self.ovl.execute_run_sensitivities(funcs, stab_derivs=stab_derivs )
+        
+        func = "CL"
+        # sens_adjoint[func][dv]
+        sens_direct = self.ovl.execute_run_sensitivities_direct(geom_dvs=geom_dvs, con_dvs=con_dvs, ref_dvs=ref_dvs, param_dvs=param_dvs,  add_stab_derivs=True)
+        import pprint
+        pprint.pprint(sens_adjoint)
+        pprint.pprint(sens_direct)
+        
+        for func in funcs:
+            for dv in geom_dvs:
+                print(f"Adjoint d{func}/d{[dv[0]]} {[dv[1]]} {sens_adjoint[func][dv[0]][dv[1]]}")
+                print(f"Direct  d{func}/d{[dv[0]]} {[dv[1]]} {sens_direct[func][dv[0]][dv[1]]}")
+            
+            for dv in con_dvs + ref_dvs + param_dvs:
+                print(f"Adjoint d{func}/d{dv} {sens_adjoint[func][dv]}")
+                print(f"Direct  d{func}/d{dv} {sens_direct[func][dv]}")
+            
 
 if __name__ == "__main__":
     unittest.main()