diff --git a/aviary/docs/examples/additional_flight_phases.ipynb b/aviary/docs/examples/additional_flight_phases.ipynb
index 293c11858..908f19636 100644
--- a/aviary/docs/examples/additional_flight_phases.ipynb
+++ b/aviary/docs/examples/additional_flight_phases.ipynb
@@ -72,6 +72,7 @@
     "            'duration_bounds': ((25.5, 76.5), 'min'),\n",
     "        },\n",
     "        'initial_guesses': {'time': ([0, 51], 'min')},\n",
+    "        \"initial_guesses\": {\"times\": ([0, 51], \"min\")},\n",
     "    },\n",
     "    'cruise_1': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -96,6 +97,9 @@
     "            'duration_bounds': ((23.5, 70.5), 'min'),\n",
     "        },\n",
     "        'initial_guesses': {'time': ([51, 47], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([51, 47], \"min\")},\n",
+
     "    },\n",
     "    'climb_2': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -119,7 +123,11 @@
     "            'initial_bounds': ((49.0, 147.0), 'min'),\n",
     "            'duration_bounds': ((5.0, 15.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([98, 10], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([98, 10], \"min\")},\n",
+
     "    },\n",
     "    'cruise_2': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -143,7 +151,11 @@
     "            'initial_bounds': ((54.0, 162.0), 'min'),\n",
     "            'duration_bounds': ((24.0, 72.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([108, 48], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([108, 48], \"min\")},\n",
+
     "    },\n",
     "    'climb_3': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -167,7 +179,11 @@
     "            'initial_bounds': ((78.0, 234.0), 'min'),\n",
     "            'duration_bounds': ((7.0, 21.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([156, 14], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([156, 14], \"min\")},\n",
+
     "    },\n",
     "    'climb_4': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -191,7 +207,11 @@
     "            'initial_bounds': ((85.0, 255.0), 'min'),\n",
     "            'duration_bounds': ((43.0, 129.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([170, 86], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([170, 86], \"min\")},\n",
+
     "    },\n",
     "    'descent_1': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -215,7 +235,11 @@
     "            'initial_bounds': ((128.0, 384.0), 'min'),\n",
     "            'duration_bounds': ((41.0, 123.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([256, 82], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([256, 82], \"min\")},\n",
+
     "    },\n",
     "    'post_mission': {\n",
     "        'include_landing': False,\n",
@@ -298,11 +322,18 @@
     "\n",
     "Playing around with a model and seeing how different settings affect the optimization and resulting aircraft design is always an enlightening experience."
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "base",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -316,9 +347,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.12.2"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }
diff --git a/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb b/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb
index 92c1d34be..daebd7818 100644
--- a/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb
+++ b/aviary/docs/examples/coupled_aircraft_mission_optimization.ipynb
@@ -68,7 +68,11 @@
     "            'initial_bounds': ((0.0, 0.0), 'min'),\n",
     "            'duration_bounds': ((35.0, 105.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([0, 70], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([0, 70], \"min\")},\n",
+
     "    },\n",
     "    'cruise': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -92,7 +96,11 @@
     "            'initial_bounds': ((35.0, 105.0), 'min'),\n",
     "            'duration_bounds': ((91.5, 274.5), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([70, 183], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([70, 183], \"min\")},\n",
+
     "    },\n",
     "    'descent_1': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -116,7 +124,11 @@
     "            'initial_bounds': ((126.5, 379.5), 'min'),\n",
     "            'duration_bounds': ((25.0, 75.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([253, 50], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([253, 50], \"min\")},\n",
+
     "    },\n",
     "    'post_mission': {\n",
     "        'include_landing': False,\n",
@@ -588,7 +600,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.12.2"
   }
  },
  "nbformat": 4,
diff --git a/aviary/docs/examples/more_advanced_example.ipynb b/aviary/docs/examples/more_advanced_example.ipynb
index 1837a19fb..d366d42ed 100644
--- a/aviary/docs/examples/more_advanced_example.ipynb
+++ b/aviary/docs/examples/more_advanced_example.ipynb
@@ -67,7 +67,11 @@
     "            'initial_bounds': ((0.0, 0.0), 'min'),\n",
     "            'duration_bounds': ((27.0, 81.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([0, 54], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([0, 54], \"min\")},\n",
+
     "    },\n",
     "    'cruise': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -91,7 +95,11 @@
     "            'initial_bounds': ((27.0, 81.0), 'min'),\n",
     "            'duration_bounds': ((85.5, 256.5), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([54, 171], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([54, 171], \"min\")},\n",
+
     "    },\n",
     "    'descent_1': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -115,7 +123,11 @@
     "            'initial_bounds': ((112.5, 337.5), 'min'),\n",
     "            'duration_bounds': ((26.5, 79.5), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([225, 53], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([225, 53], \"min\")},\n",
+
     "    },\n",
     "    'post_mission': {\n",
     "        'include_landing': False,\n",
@@ -468,7 +480,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.8.19"
   }
  },
  "nbformat": 4,
diff --git a/aviary/docs/examples/simple_mission_example.ipynb b/aviary/docs/examples/simple_mission_example.ipynb
index 7fb77f491..bdfe2b3c6 100644
--- a/aviary/docs/examples/simple_mission_example.ipynb
+++ b/aviary/docs/examples/simple_mission_example.ipynb
@@ -190,7 +190,11 @@
     "            'initial_bounds': ((0.0, 0.0), 'min'),\n",
     "            'duration_bounds': ((27.0, 81.0), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([0, 54], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([0, 54], \"min\")},\n",
+
     "    },\n",
     "    'cruise': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -214,7 +218,11 @@
     "            'initial_bounds': ((27.0, 81.0), 'min'),\n",
     "            'duration_bounds': ((85.5, 256.5), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([54, 171], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([54, 171], \"min\")},\n",
+
     "    },\n",
     "    'descent_1': {\n",
     "        'subsystem_options': {'core_aerodynamics': {'method': 'computed'}},\n",
@@ -238,7 +246,11 @@
     "            'initial_bounds': ((112.5, 337.5), 'min'),\n",
     "            'duration_bounds': ((26.5, 79.5), 'min'),\n",
     "        },\n",
+
     "        'initial_guesses': {'time': ([225, 53], 'min')},\n",
+
+    "        \"initial_guesses\": {\"times\": ([225, 53], \"min\")},\n",
+
     "    },\n",
     "    'post_mission': {\n",
     "        'include_landing': False,\n",
@@ -418,7 +430,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.8.19"
   }
  },
  "nbformat": 4,
diff --git a/aviary/mission/sixdof/force_component_calc.py b/aviary/mission/sixdof/force_component_calc.py
new file mode 100644
index 000000000..9fe759467
--- /dev/null
+++ b/aviary/mission/sixdof/force_component_calc.py
@@ -0,0 +1,490 @@
+import numpy as np
+import openmdao.api as om
+from aviary.variable_info.variables import Dynamic
+from aviary.utils.functions import add_aviary_input
+
+
+class ForceComponentResolver(om.ExplicitComponent):
+    """
+    This class will resolve forces (thrust, drag, lift, etc.) into their 
+    respective x,y,z components for the 6 DOF equations of motion. 
+    
+    This class assumes that the total force is given and needs to be resolved 
+    into the separate components.
+
+    Assumptions:
+        - Thrust is entirely in -z direction (T = (0,0,-T_z)^T) w.r.t. body CS
+        - Assuming F_i is in body CS, and D, S, and L are in wind CS. Wind -> body rotation matrix
+          was applied for coordinate transformations
+        - Thrust is initially in NED CS. So, two rotations (NED -> wind and wind -> body) are applied
+
+    """
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        # inputs
+
+        self.add_input(
+            'u',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="axial velocity"
+        )
+
+        self.add_input(
+            'v',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="lateral velocity"
+        )
+
+        self.add_input(
+            'w',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="vertical velocity"
+        )
+
+        self.add_input(
+            'drag',
+            val=np.zeros(nn),
+            units='N',
+            desc="Drag vector (unresolved)"
+        )
+
+        self.add_input(
+            'thrust',
+            val=np.zeros(nn),
+            units='N',
+            desc="Thrust vector (unresolved)"
+        )
+
+        self.add_input(
+            'lift',
+            val=np.zeros(nn),
+            units='N',
+            desc="Lift vector (unresolved)"
+        )
+
+        self.add_input(
+            'side',
+            val=np.zeros(nn),
+            units='N',
+            desc="Side vector (unresolved)"
+        )
+
+        self.add_input(
+            'heading_angle',
+            val=np.zeros(nn),
+            units='rad',
+            desc='Heading angle in body'
+        )
+
+        add_aviary_input(self,
+                         Dynamic.Mission.FLIGHT_PATH_ANGLE,
+                         units='rad',
+                         desc="Flight path angle in body")
+        
+        self.add_input(
+            'heading_angle_NED',
+            val=np.zeros(nn),
+            units='rad',
+            desc="Thrust heading angle in NED"
+        )
+
+        self.add_input(
+            'fpa_NED',
+            val=np.zeros(nn),
+            units='rad',
+            desc="Thrust flight path angle in NED"
+        )
+
+        # self.add_input(
+        #     'true_air_speed',
+        #     val=np.zeros(nn),
+        #     units='m/s',
+        #     desc="True air speed"
+        # ) # This is an aviary variable
+
+        # outputs
+
+        self.add_output(
+            'Fx',
+            val=np.zeros(nn),
+            units='N',
+            desc="x-comp of final force"
+        )
+
+        self.add_output(
+            'Fy',
+            val=np.zeros(nn),
+            units='N',
+            desc="y-comp of final force"
+        )
+
+        self.add_output(
+            'Fz',
+            val=np.zeros(nn),
+            units='N',
+            desc="z-comp of final force"
+        )
+
+        ar = np.arange(nn)
+
+        self.declare_partials(of='Fx', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='drag', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='lift', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='side', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='thrust', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='heading_angle', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt=Dynamic.Mission.FLIGHT_PATH_ANGLE, rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='heading_angle_NED', rows=ar, cols=ar)
+        self.declare_partials(of='Fx', wrt='fpa_NED', rows=ar, cols=ar)
+
+        self.declare_partials(of='Fy', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='drag', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='lift', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='side', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='thrust', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='heading_angle', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt=Dynamic.Mission.FLIGHT_PATH_ANGLE, rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='heading_angle_NED', rows=ar, cols=ar)
+        self.declare_partials(of='Fy', wrt='fpa_NED', rows=ar, cols=ar)
+
+        self.declare_partials(of='Fz', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='drag', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='lift', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='side', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='thrust', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='heading_angle', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt=Dynamic.Mission.FLIGHT_PATH_ANGLE, rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='heading_angle_NED', rows=ar, cols=ar)
+        self.declare_partials(of='Fz', wrt='fpa_NED', rows=ar, cols=ar)
+
+    def compute(self, inputs, outputs):
+
+        u = inputs['u']
+        v = inputs['v']
+        w = inputs['w']
+        D = inputs['drag']
+        T = inputs['thrust']
+        L = inputs['lift']
+        S = inputs['side'] # side force -- assume 0 for now
+        chi = inputs['heading_angle']
+        gamma = inputs[Dynamic.Mission.FLIGHT_PATH_ANGLE]
+        chi_T = inputs['heading_angle_NED']
+        gamma_T = inputs['fpa_NED']
+
+        nn = self.options['num_nodes']
+
+        # true air speed
+
+        V = np.sqrt(u**2 + v**2 + w**2)
+
+        # angle of attack
+
+        # divide by zero checks
+        if np.any(u == 0):
+            u[u == 0] = 1e-4
+            alpha = np.arctan(w / u)
+        else:
+            alpha = np.arctan(w / u)
+
+        # side slip angle
+
+        # divide by zero checks
+        if ((np.any(u != 0) or np.any(w != 0))) :
+            beta = np.arctan(v / np.sqrt(u**2 + w**2))
+        else:
+            u[u == 0] = 1.0e-4
+            beta = np.arctan(v / np.sqrt(u**2 + w**2))
+        
+        
+        
+
+        # some trig needed
+
+        cos_a = np.cos(alpha)
+        cos_b = np.cos(beta)
+        sin_a = np.sin(alpha)
+        sin_b = np.sin(beta)
+        cos_g = np.cos(gamma)
+        sin_g = np.sin(gamma)
+        cos_c = np.cos(chi)
+        sin_c = np.sin(chi)
+
+        # Thrust direction in NED -- as \hat{t}_n
+
+        t_hat_n = [
+            np.cos(gamma_T) * np.cos(chi_T),
+            np.cos(gamma_T) * np.sin(chi_T),
+            -np.sin(gamma_T)
+        ]
+
+        t_hat_n = np.array(t_hat_n)
+        t_hat_n = t_hat_n.reshape((3, 1))
+
+        # C_{b-<n}
+
+        Cbn = [
+            -sin_b * sin_c - cos_b * cos_c * np.cos(alpha - gamma),
+             sin_b * cos_c - cos_b * sin_c * np.cos(alpha - gamma),
+             cos_b * np.sin(alpha - gamma),
+            cos_b * sin_c - sin_b * cos_c * np.cos(alpha - gamma),
+             -cos_b * cos_c - sin_b * sin_c * np.cos(alpha - gamma),
+             sin_b * np.sin(alpha - gamma),
+            cos_c * np.sin(alpha - gamma),
+             sin_c * np.sin(alpha - gamma),
+             np.cos(alpha - gamma)
+        ]
+
+        Cbn = np.array(Cbn)
+        Cbn = Cbn.reshape((3, 3))
+
+        
+
+        print("Cbn shape: ", np.shape(Cbn))
+        print("t_hat_n shape: ", np.shape(t_hat_n))
+        print("Cbn = ", Cbn)
+        print("t_hat_n = ", t_hat_n)
+        
+
+        # Thrust direction in body
+        t_hat_b = Cbn @ t_hat_n
+
+        # Thrust force in body
+        Fb_thrust = T * t_hat_b
+
+        # Thrust in (x,y,z)
+        F_T_x = Fb_thrust[0]
+        F_T_y = Fb_thrust[1]
+        F_T_z = Fb_thrust[2]
+        
+
+        outputs['Fx'] = -(cos_a * cos_b * D - cos_a * sin_b * S - sin_a * L) + F_T_x
+        outputs['Fy'] = -(sin_b * D + cos_b * S) + F_T_y
+        outputs['Fz'] = -(sin_a * cos_b * D + sin_a * sin_b * S + cos_a * L) + F_T_z
+
+    def compute_partials(self, inputs, J):
+        u = inputs['u']
+        v = inputs['v']
+        w = inputs['w']
+        D = inputs['drag']
+        T = inputs['thrust']
+        L = inputs['lift']
+        S = inputs['side'] # side force -- assume 0 for now
+        gamma = inputs[Dynamic.Mission.FLIGHT_PATH_ANGLE]
+        chi = inputs['heading_angle']
+        chi_T = inputs['heading_angle_NED']
+        gamma_T = inputs['fpa_NED']
+
+        V = np.sqrt(u**2 + v**2 + w**2)
+
+        # divide by zero checks
+        if u == 0:
+            u = 1e-4
+            alpha = np.arctan(w / u)
+        else:
+            alpha = np.arctan(w / u)
+
+        # side slip angle
+
+        # divide by zero checks
+        if ((u != 0 or w != 0)) :
+            beta = np.arctan(v / np.sqrt(u**2 + w**2))
+        else:
+            u = 1.0e-4
+            beta = np.arctan(v / np.sqrt(u**2 + w**2))
+        
+        t_hat_n = [
+            np.cos(gamma_T) * np.cos(chi_T),
+            np.cos(gamma_T) * np.sin(chi_T),
+            -np.sin(gamma_T)
+        ]
+
+        t_hat_n = np.array(t_hat_n)
+        t_hat_n = t_hat_n.reshape((3, 1))
+
+        cos_a, sin_a = np.cos(alpha), np.sin(alpha)
+        cos_b, sin_b = np.cos(beta),  np.sin(beta)
+        cos_g, sin_g = np.cos(gamma), np.sin(gamma)
+        cos_c, sin_c = np.cos(chi),   np.sin(chi)
+
+        # C_{b-<n}
+
+        Cbn = [
+            -sin_b * sin_c - cos_b * cos_c * np.cos(alpha - gamma),
+             sin_b * cos_c - cos_b * sin_c * np.cos(alpha - gamma),
+             cos_b * np.sin(alpha - gamma),
+            cos_b * sin_c - sin_b * cos_c * np.cos(alpha - gamma),
+             -cos_b * cos_c - sin_b * sin_c * np.cos(alpha - gamma),
+             sin_b * np.sin(alpha - gamma),
+            cos_c * np.sin(alpha - gamma),
+             sin_c * np.sin(alpha - gamma),
+             np.cos(alpha - gamma)
+        ]
+
+        Cbn = np.array(Cbn)
+        Cbn = Cbn.reshape((3, 3))
+
+        # Derivatives of t_n
+        dt_n_dchi = [
+            -np.cos(gamma_T) * np.sin(chi_T),
+            np.cos(gamma_T) * np.cos(chi_T),
+            0.0
+        ]
+
+        dt_n_dchi = np.array(dt_n_dchi, dtype="object")
+        dt_n_dchi = dt_n_dchi.reshape((3, 1))
+
+        dt_n_dgamma = [
+            -np.sin(gamma_T) * np.cos(chi_T),
+            -np.sin(gamma_T) * np.sin(chi_T),
+            -np.cos(gamma_T)
+        ]
+
+        dt_n_dgamma = np.array(dt_n_dgamma, dtype="object")
+        dt_n_dgamma = dt_n_dgamma.reshape((3, 1))
+
+        # Rotation to body + derivatives
+        t_b = Cbn @ t_hat_n
+        dt_b_dchi = Cbn @ dt_n_dchi
+        dt_b_dgamma = Cbn @ dt_n_dgamma
+
+        # note: d/dx arctan(a/x) = -a / (x^2 + a^2)
+        # note: d/dx arctan(a / sqrt(b^2 + x^2)) = - ax / ((b^2 + x^2 + a^2) * sqrt(b^2 + x^2))
+        # note: d/dx arctan(x / sqrt(a^2 + b^2)) = sqrt(a^2 + b^2) / (a^2 + b^2 + x^2)
+        # note: d/dx arctan(x/a) = a / (a^2 + x^2)
+
+        J['Fx', 'u'] = np.cos(alpha) * np.sin(beta) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * D + \
+                         np.cos(beta) * np.sin(alpha) * (-w / (w**2 + u**2)) * D + \
+                         (np.cos(alpha) * np.cos(beta) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * S - np.sin(beta) * np.sin(alpha) * (-w / (w**2 + u**2)) * S) + \
+                         (np.cos(alpha) * (-w / (w**2 + u**2)) * L) + T * ((-np.cos(gamma_T) * np.cos(chi_T) * cos_b * sin_c * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) + 
+                                                                           np.cos(gamma_T) * np.cos(chi_T) * sin_b * cos_c * np.cos(alpha - gamma) * 
+                                                                           ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) + np.cos(gamma_T) * np.cos(chi_T) * cos_b * cos_c * 
+                                                                           np.sin(alpha - gamma) * (-w / (w**2 + u**2))) + (
+                                                                               np.cos(gamma_T) * np.sin(chi_T) * cos_b * cos_c * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) + 
+                                                                                np.cos(gamma_T) * np.sin(chi_T) * sin_b * sin_c * np.cos(alpha - gamma) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) + 
+                                                                           np.cos(gamma_T) * np.sin(chi_T) * cos_b * sin_c * np.sin(alpha -gamma) * (-w / (w**2 + u**2))) + 
+                                                                           (np.sin(gamma_T) * sin_b * np.sin(alpha - gamma) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) - 
+                                                                            np.sin(gamma_T) * cos_b * np.cos(alpha - gamma) * (-w / (w**2 + u**2))))
+        J['Fx', 'v'] = np.cos(alpha) * np.sin(beta) * (np.sqrt(w**2 + u**2) / V**2) * D + np.cos(alpha) * np.cos(beta) * (np.sqrt(w**2 + u**2) / V**2) * S + T * (
+            (-np.cos(gamma_T) * np.cos(chi_T) * cos_b * sin_c * (np.sqrt(w**2 + u**2) / V**2) + np.cos(gamma_T) * np.cos(chi_T) * sin_b * cos_c * np.cos(alpha - gamma) * (np.sqrt(w**2 + u**2) / V**2)) + 
+             (np.cos(gamma_T) * np.sin(chi_T) * cos_b * cos_c * (np.sqrt(w**2 + u**2) / V**2) + np.cos(gamma_T) * np.sin(chi_T) * sin_b * sin_c * np.cos(alpha - gamma) * (np.sqrt(w**2 + u**2) / V**2)) +
+             (np.sin(gamma_T) * sin_b * np.sin(alpha - gamma) * (np.sqrt(w**2 + u**2) / V**2))
+        )
+        J['Fx', 'w'] = np.cos(alpha) * np.sin(beta) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * D + np.sin(alpha) * np.cos(beta) * (u / (w**2 + u**2)) * D + \
+                       np.cos(alpha) * np.cos(beta) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * S - np.sin(alpha) * np.sin(beta) * (u / (w**2 + u**2)) * S + \
+                       np.cos(alpha) * (u / (w**2 + u**2)) * L + T * (
+                           (-np.cos(gamma_T) * np.cos(chi_T) * cos_b * sin_c * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) + np.cos(gamma_T) * np.cos(chi_T) * cos_b * cos_c * np.sin(alpha - gamma) * (u / (w**2 + u**2)) +
+                            np.cos(gamma_T) * np.cos(chi_T) * sin_b * cos_c * np.cos(alpha - gamma) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2))))) + 
+                           (np.cos(gamma_T) * np.sin(chi_T) * cos_b * cos_c * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) + np.cos(gamma_T) * np.sin(chi_T) * cos_b * sin_c * np.sin(alpha - gamma) * (u / (w**2 + u**2)) + 
+                            np.cos(gamma_T) * np.sin(chi_T) * sin_b * sin_c * np.cos(alpha - gamma) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2))))) + 
+                           (-np.sin(gamma_T) * cos_b * np.cos(alpha - gamma) * (u / (w**2 + u**2)) + np.sin(gamma_T) * sin_b * np.sin(alpha - gamma) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))))
+                       )
+        J['Fx', 'drag'] = -np.cos(alpha) * np.cos(beta)
+        J['Fx', 'lift'] = np.sin(alpha)
+        J['Fx', 'side'] = np.cos(alpha) * np.sin(beta)
+        J['Fx', 'thrust'] = t_b[0]
+        J['Fx', Dynamic.Mission.FLIGHT_PATH_ANGLE] = T * ((-cos_b * cos_c * np.sin(alpha - gamma)) * t_hat_n[0] + (-cos_b * sin_c * np.sin(alpha - gamma)) * t_hat_n[1] + 
+                                                          (-cos_b * np.cos(alpha - gamma)) * t_hat_n[2])
+        J['Fx', 'heading_angle'] = T * ((-sin_b * cos_c + cos_b * sin_c * np.cos(alpha - gamma)) * t_hat_n[0] + 
+                                        (-sin_b * sin_c - cos_b * cos_c * np.cos(alpha - gamma)) * t_hat_n[1])
+        J['Fx', 'heading_angle_NED'] = T * dt_b_dchi[0]
+        J['Fx', 'fpa_NED'] = T * dt_b_dgamma[0]
+
+        J['Fy', 'u'] = -np.cos(beta) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * D + np.sin(beta) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * S - \
+                        T * np.cos(gamma_T) * np.cos(chi_T) * sin_b * sin_c * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) - \
+                        np.cos(gamma_T) * np.cos(chi_T) * np.cos(beta) * cos_c * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * T * np.cos(alpha - gamma) + \
+                        np.sin(alpha - gamma) * cos_c * (-w / (w**2 + u**2)) * T * np.sin(beta) * np.cos(gamma_T) * np.cos(chi_T) + \
+                        sin_b * cos_c * np.cos(gamma_T) * np.sin(chi_T) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * T - \
+                        np.cos(gamma_T) * np.sin(chi_T) * cos_b * sin_c * np.cos(alpha - gamma) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * T + \
+                        np.cos(gamma_T) * np.sin(chi_T) * sin_b * sin_c * np.sin(alpha - gamma) * (-w / (w**2 + u**2)) * T - \
+                        np.sin(gamma_T) * cos_b * np.sin(alpha - gamma) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * T - \
+                        np.sin(gamma_T) * sin_b * np.cos(alpha - gamma) * (-w / (w**2 + u**2)) * T
+        J['Fy', 'v'] = -np.cos(beta) * (np.sqrt(w**2 + u**2) / V**2) * D + np.sin(beta) * (np.sqrt(w**2 + u**2) / V**2) * S - \
+                        np.cos(gamma_T) * np.cos(chi_T) * sin_b * sin_c * (np.sqrt(w**2 + u**2) / V**2) * T - \
+                        np.cos(gamma_T) * np.cos(chi_T) * cos_b * cos_c * np.cos(alpha - gamma) * (np.sqrt(w**2 + u**2) / V**2) * T + \
+                        np.cos(gamma_T) * np.sin(chi_T) * sin_b * cos_c * (np.sqrt(w**2 + u**2) / V**2) * T - \
+                        np.cos(gamma_T) * np.sin(chi_T) * np.cos(beta) * sin_c * (np.sqrt(w**2 + u**2) / V**2) * T * np.cos(alpha - gamma) - \
+                        np.sin(gamma_T) * cos_b * np.sin(alpha - gamma) * (np.sqrt(w**2 + u**2) / V**2) * T
+        J['Fy', 'w'] = -np.cos(beta) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * D + np.sin(beta) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * S - \
+                        np.cos(gamma_T) * np.cos(chi_T) * sin_b * sin_c * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * T - \
+                        np.cos(gamma_T) * np.cos(chi_T) * cos_b * cos_c * np.cos(alpha - gamma) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * T + \
+                        np.cos(gamma_T) * np.cos(chi_T) * sin_b * cos_c * np.sin(alpha - gamma) * (u / (w**2 + u**2)) * T  + \
+                        np.cos(gamma_T) * np.sin(chi_T) * sin_b * cos_c * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * T - \
+                        np.cos(gamma_T) * np.sin(chi_T) * cos_b * sin_c * np.cos(alpha - gamma) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * T + \
+                        np.cos(gamma_T) * np.sin(chi_T) * sin_b * sin_c * np.sin(alpha - gamma) * (u / (w**2 + u**2)) * T - \
+                        np.sin(gamma_T) * np.cos(beta) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * T * np.sin(alpha - gamma) - \
+                        np.sin(gamma_T) * np.cos(alpha - gamma) * (u / (w**2 + u**2)) * T * np.sin(beta)
+        J['Fy', 'drag'] = -np.sin(beta)
+        J['Fy', 'side'] = -np.cos(beta)
+        J['Fy', 'thrust'] = t_b[1]
+        J['Fy', Dynamic.Mission.FLIGHT_PATH_ANGLE] = T * ((-sin_b * cos_c * np.sin(alpha - gamma)) * t_hat_n[0] + (-sin_b * sin_c * np.sin(alpha - gamma)) * t_hat_n[1] + 
+                                                          (-sin_b * np.cos(alpha - gamma)) * t_hat_n[2])
+        J['Fy', 'heading_angle'] = T * ((cos_b * cos_c + sin_b * sin_c * np.cos(alpha - gamma)) * t_hat_n[0] + (cos_b * sin_c - sin_b * cos_c * np.cos(alpha - gamma)) * t_hat_n[1])
+        J['Fy', 'heading_angle_NED'] = T * dt_b_dchi[1]
+        J['Fy', 'fpa_NED'] = T * dt_b_dgamma[1]
+
+        J['Fz', 'u'] = np.sin(alpha) * np.sin(beta) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * D - np.cos(alpha) * np.cos(beta) * (-w / (w**2 + u**2)) * D - \
+                       np.sin(alpha) * np.cos(beta) * ((-v * u) / ((V**2 * np.sqrt(w**2 + u**2)))) * S - np.cos(alpha) * np.sin(beta) * (-w / (w**2 + u**2)) * S + \
+                       np.sin(alpha) * (-w / (w**2 + u**2)) * L + np.cos(gamma_T) * np.cos(chi_T) * cos_c * np.cos(alpha - gamma) * (-w / (w**2 + u**2)) * T + \
+                       np.cos(gamma_T) * np.sin(chi_T) * sin_c * np.cos(alpha - gamma) * (-w / (w**2 + u**2)) * T + \
+                       np.sin(gamma_T) * np.sin(alpha - gamma) * (-w / (w**2 + u**2)) * T
+        J['Fz', 'v'] = np.sin(alpha) * np.sin(beta) * (np.sqrt(w**2 + u**2) / V**2) * D - np.sin(alpha) * np.cos(beta) * (np.sqrt(w**2 + u**2) / V**2) * S 
+        J['Fz', 'w'] = np.sin(alpha) * np.sin(beta) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * D - np.cos(alpha) * np.cos(beta) * (u / (w**2 + u**2)) * D - \
+                       np.sin(alpha) * np.cos(beta) * ((-w * v) / ((V**2 * np.sqrt(w**2 + u**2)))) * S - np.cos(alpha) * np.sin(beta) * (u / (w**2 + u**2)) * S + \
+                       np.sin(alpha) * (u / (w**2 + u**2)) * L + np.cos(gamma_T) * np.cos(chi_T) * cos_c * np.cos(alpha - gamma) * (u / (w**2 + u**2)) * T + \
+                       np.cos(gamma_T) * np.sin(chi_T) * sin_c * np.cos(alpha - gamma) * (u / (w**2 + u**2)) * T + \
+                       np.sin(gamma_T) * np.sin(alpha - gamma) * (u / (w**2 + u**2)) * T
+        J['Fz', 'drag'] = -np.sin(alpha) * np.cos(beta) 
+        J['Fz', 'lift'] = -np.cos(alpha)
+        J['Fz', 'side'] = -np.sin(alpha) * np.sin(beta)
+        J['Fz', 'thrust'] = t_b[2]
+        J['Fz', Dynamic.Mission.FLIGHT_PATH_ANGLE] = T * ((-cos_c * np.cos(alpha - gamma)) * t_hat_n[0] + (-sin_c * np.cos(alpha - gamma)) * t_hat_n[1] + 
+                                                          np.sin(alpha - gamma) * t_hat_n[2])
+        J['Fz', 'heading_angle'] = T * ((-sin_c * np.sin(alpha - gamma)) * t_hat_n[0] + (cos_c * np.sin(alpha - gamma)) * t_hat_n[1])
+        J['Fz', 'heading_angle_NED'] = T * dt_b_dchi[2]
+        J['Fz', 'fpa_NED'] = T * dt_b_dgamma[2]
+
+if __name__ == "__main__":
+    p = om.Problem()
+    p.model = om.Group()
+    des_vars = p.model.add_subsystem('des_vars', om.IndepVarComp(), promotes=['*'])
+
+    des_vars.add_output('u', 0.5, units='m/s')
+    des_vars.add_output('v', 0.6, units='m/s')
+    des_vars.add_output('w', 0.7, units='m/s')
+    des_vars.add_output('drag', 50, units='N')
+    des_vars.add_output('thrust', 50, units='N')
+    des_vars.add_output('lift', 60, units='N')
+    des_vars.add_output('side', 70, units='N')
+    des_vars.add_output('heading_angle', units='rad')
+    des_vars.add_output(Dynamic.Mission.FLIGHT_PATH_ANGLE, units='rad')
+    des_vars.add_output('heading_angle_NED', units='rad')
+    des_vars.add_output('fpa_NED', units='rad')
+
+    p.model.add_subsystem('ForceComponentResolver', ForceComponentResolver(num_nodes=1), promotes=['*'])
+
+    p.setup(check=False, force_alloc_complex=True)
+
+    p.run_model()
+
+    p.check_partials(compact_print=True, show_only_incorrect=True, method='cs')
+
+
+
+
+
+
+
diff --git a/aviary/mission/sixdof/six_dof_EOM.py b/aviary/mission/sixdof/six_dof_EOM.py
new file mode 100644
index 000000000..55e8adc25
--- /dev/null
+++ b/aviary/mission/sixdof/six_dof_EOM.py
@@ -0,0 +1,738 @@
+import numpy as np
+import openmdao.api as om
+
+class SixDOF_EOM(om.ExplicitComponent):
+    """
+    Six DOF EOM component, with particular emphasis for rotorcraft. 
+    ASSUMPTIONS:
+        - Assume Flat Earth model (particularly for rotorcraft)
+        - Earth is the internal f.o.r.
+        - Gravity is constant and normal to the tangent plane (Earth's surface) -- g = (0 0 G)^T
+        - (aircraft) mass is constant
+        - aircraft is a rigid body
+        - Symmetry in the xz plane (for moment of inertia matrix -- J_xy = J_yx = J_zy = J_yz = 0)
+
+    """
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+    
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        self.add_input(
+            'mass',
+            val=np.zeros(nn),
+            units='kg',
+            desc="mass -- assume constant"
+        )
+
+        self.add_input(
+            'u',
+            val=np.zeros(nn),
+            units='m/s', # meters per second
+            desc="axial velocity of CM wrt inertial CS resolved in aircraft body fixed CS"
+        )
+
+        self.add_input(
+            'v',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="lateral velocity of CM wrt inertial CS resolved in aircraft body fixed CS"
+        )
+
+        self.add_input(
+            'w',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="vertical velocity of CM wrt inertial CS resolved in aircraft body fixed CS"
+        )
+
+        self.add_input(
+            'roll_ang_vel',
+            val=np.zeros(nn),
+            units='rad/s', # radians per second
+            desc="roll angular velocity of body fixed CS wrt intertial CS"
+        )
+
+        self.add_input(
+            'pitch_ang_vel',
+            val=np.zeros(nn),
+            units='rad/s',
+            desc="pitch angular velocity of body fixed CS wrt intertial CS"
+        )
+
+        self.add_input(
+            'yaw_ang_vel',
+            val=np.zeros(nn),
+            units='rad/s',
+            desc="yaw angular velocity of body fixed CS wrt intertial CS"
+        )
+
+        self.add_input(
+            'roll',
+            val=np.zeros(nn),
+            units='rad', # radians
+            desc="roll angle"
+        )
+
+        self.add_input(
+            'pitch',
+            val=np.zeros(nn),
+            units='rad',
+            desc="pitch angle"
+        )
+
+        self.add_input(
+            'yaw',
+            val=np.zeros(nn),
+            units='rad',
+            desc="yaw angle"
+        )
+
+        self.add_input(
+            'x',
+            val=np.zeros(nn),
+            units='m',
+            desc="x-axis position of aircraft resolved in North-East-Down (NED) CS"
+        )
+
+        self.add_input(
+            'y',
+            val=np.zeros(nn),
+            units='m',
+            desc="y-axis position of aircraft resolved in NED CS"
+        )
+
+        self.add_input(
+            'z',
+            val=np.zeros(nn),
+            units='m',
+            desc="z-axis position of aircraft resolved in NED CS"
+        )
+
+        self.add_input(
+            'g',
+            val=9.81,
+            units='m/s**2',
+            desc="acceleration due to gravity"
+        )
+
+        self.add_input(
+            'Fx',
+            val=np.zeros(nn),
+            units='N',
+            desc="external forces in the x direciton"
+        )
+
+        self.add_input(
+            'Fy',
+            val=np.zeros(nn),
+            units='N',
+            desc="external forces in the y direction"
+        )
+
+        self.add_input(
+            'Fz',
+            val=np.zeros(nn),
+            units='N',
+            desc="external forces in the z direction"
+        )
+
+        self.add_input(
+            'lx',
+            val=np.zeros(nn),
+            units='kg*m**2/s**2', # kg times m^2 / s^2
+            desc="external moments in the x direction"
+        )
+
+        self.add_input(
+            'ly',
+            val=np.zeros(nn),
+            units='kg*m**2/s**2',
+            desc="external moments in the y direction"
+        )
+
+        self.add_input(
+            'lz',
+            val=np.zeros(nn),
+            units='kg*m**2/s**2',
+            desc="external moments in the z direction"
+        )
+
+        # Below are the necessary components for the moment of inertia matrix (J)
+        # Only xx, yy, zz, and xz are needed (xy and yz are 0 with assumptions)
+        # For now, these are separated.
+        # TODO: Rewrite J and EOM in matrix form
+
+        self.add_input(
+            'J_xz',
+            val=np.zeros(1),
+            units='kg*m**2',
+            desc="x-z (top right and bottom left corner of 3x3 matrix, assuming symmetry) " \
+            "component"
+        )
+
+        self.add_input(
+            'J_xx',
+            val=np.zeros(1),
+            units='kg*m**2',
+            desc="first diag component"
+        )
+
+        self.add_input(
+            'J_yy',
+            val=np.zeros(1),
+            units='kg*m**2',
+            desc="second diag component"
+        )
+
+        self.add_input(
+            'J_zz',
+            val=np.zeros(1),
+            units='kg*m**2',
+            desc="third diag component"
+        )
+
+        # Outputs
+
+        self.add_output(
+            'dx_accel',
+            val=np.zeros(nn),
+            units='m/s**2', # meters per seconds squared
+            desc="x-axis (roll-axis) velocity equation, " \
+            "state: u",
+            tags=['dymos.state_rate_source:u', 'dymos.state_units:m/s']
+        )
+
+        self.add_output(
+            'dy_accel',
+            val=np.zeros(nn),
+            units='m/s**2',
+            desc="y-axis (pitch axis) velocity equation, " \
+            "state: v",
+            tags=['dymos.state_rate_source:v', 'dymos.state_units:m/s']
+        )
+
+        self.add_output(
+            'dz_accel',
+            val=np.zeros(nn),
+            units='m/s**2',
+            desc="z-axis (yaw axis) velocity equation, " \
+            "state: w",
+            tags=['dymos.state_rate_source:w', 'dymos.state_units:m/s']
+        )
+
+        self.add_output(
+            'roll_accel',
+            val=np.zeros(nn),
+            units='rad/s**2', # radians per second squared
+            desc="roll equation, " \
+            "state: roll_ang_vel",
+            tags=['dymos.state_rate_source:roll_ang_vel', 'dymos.state_units:rad/s']
+        )
+
+        self.add_output(
+            'pitch_accel',
+            val=np.zeros(nn),
+            units='rad/s**2',
+            desc="pitch equation, " \
+            "state: pitch_ang_vel",
+            tags=['dymos.state_rate_source:pitch_ang_vel', 'dymos.state_units:rad/s']
+        )
+
+        self.add_output(
+            'yaw_accel',
+            val=np.zeros(nn),
+            units='rad/s**2',
+            desc="yaw equation, " \
+            "state: yaw_ang_vel",
+            tags=['dymos.state_rate_source:yaw_ang_vel', 'dymos.state_units:rad/s']
+        )
+
+        self.add_output(
+            'roll_angle_rate_eq',
+            val=np.zeros(nn),
+            units='rad/s',
+            desc="Euler angular roll rate",
+            tags=['dymos.state_rate_source:roll', 'dymos.state_units:rad']
+        )
+
+        self.add_output(
+            'pitch_angle_rate_eq',
+            val=np.zeros(nn),
+            units='rad/s',
+            desc="Euler angular pitch rate",
+            tags=['dymos.state_rate_source:pitch', 'dymos.state_units:rad']
+        )
+
+        self.add_output(
+            'yaw_angle_rate_eq',
+            val=np.zeros(nn),
+            units='rad/s',
+            desc="Euler angular yaw rate",
+            tags=['dymos.state_rate_source:yaw', 'dymos.state_units:rad']
+        )
+
+        self.add_output(
+            'dx_dt',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="x-position derivative of aircraft COM wrt point in NED CS",
+            tags=['dymos.state_rate_source:x', 'dymos.state_units:m'] 
+        )
+
+        self.add_output(
+            'dy_dt',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="y-position derivative of aircraft COM wrt point in NED CS",
+            tags=['dymos.state_rate_source:y', 'dymos.state_units:m']
+        )
+
+        self.add_output(
+            'dz_dt',
+            val=np.zeros(nn),
+            units='m/s',
+            desc="z-position derivative of aircraft COM wrt point in NED CS",
+            tags=['dymos.state_rate_source:z', 'dymos.state_units:m']
+        )
+
+        ar = np.arange(nn)
+        self.declare_partials(of='dx_accel', wrt='mass', rows=ar, cols=ar)
+        self.declare_partials(of='dx_accel', wrt='Fx', rows=ar, cols=ar)
+        self.declare_partials(of='dx_accel', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='dx_accel', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='dx_accel', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='dx_accel', wrt='pitch_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='dx_accel', wrt='g')
+        self.declare_partials(of='dx_accel', wrt='pitch', rows=ar, cols=ar)
+
+        self.declare_partials(of='dy_accel', wrt='mass', rows=ar, cols=ar)
+        self.declare_partials(of='dy_accel', wrt='Fy', rows=ar, cols=ar)
+        self.declare_partials(of='dy_accel', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='dy_accel', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='dy_accel', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='dy_accel', wrt='roll_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='dy_accel', wrt='g')
+        self.declare_partials(of='dy_accel', wrt='roll', rows=ar, cols=ar)
+        self.declare_partials(of='dy_accel', wrt='pitch', rows=ar, cols=ar)
+
+        self.declare_partials(of='dz_accel', wrt='mass', rows=ar, cols=ar)
+        self.declare_partials(of='dz_accel', wrt='Fz', rows=ar, cols=ar)
+        self.declare_partials(of='dz_accel', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='dz_accel', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='dz_accel', wrt='roll_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='dz_accel', wrt='pitch_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='dz_accel', wrt='g')
+        self.declare_partials(of='dz_accel', wrt='roll', rows=ar, cols=ar)
+        self.declare_partials(of='dz_accel', wrt='pitch', rows=ar, cols=ar)
+
+        self.declare_partials(of='roll_accel', wrt='J_xz')
+        self.declare_partials(of='roll_accel', wrt='J_xx')
+        self.declare_partials(of='roll_accel', wrt='J_yy')
+        self.declare_partials(of='roll_accel', wrt='J_zz')
+        self.declare_partials(of='roll_accel', wrt='roll_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='roll_accel', wrt='pitch_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='roll_accel', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='roll_accel', wrt='lx', rows=ar, cols=ar)
+        self.declare_partials(of='roll_accel', wrt='lz', rows=ar, cols=ar)
+
+        self.declare_partials(of='pitch_accel', wrt='J_xz')
+        self.declare_partials(of='pitch_accel', wrt='J_xx')
+        self.declare_partials(of='pitch_accel', wrt='J_yy')
+        self.declare_partials(of='pitch_accel', wrt='J_zz')
+        self.declare_partials(of='pitch_accel', wrt='roll_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='pitch_accel', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='pitch_accel', wrt='ly', rows=ar, cols=ar)
+
+        self.declare_partials(of='yaw_accel', wrt='J_xz')
+        self.declare_partials(of='yaw_accel', wrt='J_xx')
+        self.declare_partials(of='yaw_accel', wrt='J_yy')
+        self.declare_partials(of='yaw_accel', wrt='J_zz')
+        self.declare_partials(of='yaw_accel', wrt='roll_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='yaw_accel', wrt='pitch_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='yaw_accel', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='yaw_accel', wrt='lx', rows=ar, cols=ar)
+        self.declare_partials(of='yaw_accel', wrt='lz', rows=ar, cols=ar)
+
+        self.declare_partials(of='roll_angle_rate_eq', wrt='roll_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='roll_angle_rate_eq', wrt='pitch_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='roll_angle_rate_eq', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='roll_angle_rate_eq', wrt='roll', rows=ar, cols=ar)
+        self.declare_partials(of='roll_angle_rate_eq', wrt='pitch', rows=ar, cols=ar)
+
+        self.declare_partials(of='pitch_angle_rate_eq', wrt='pitch_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='pitch_angle_rate_eq', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='pitch_angle_rate_eq', wrt='roll', rows=ar, cols=ar)
+
+        self.declare_partials(of='yaw_angle_rate_eq', wrt='pitch_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='yaw_angle_rate_eq', wrt='yaw_ang_vel', rows=ar, cols=ar)
+        self.declare_partials(of='yaw_angle_rate_eq', wrt='roll', rows=ar, cols=ar)
+        self.declare_partials(of='yaw_angle_rate_eq', wrt='pitch', rows=ar, cols=ar)
+
+        self.declare_partials(of='dx_dt', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='dx_dt', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='dx_dt', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='dx_dt', wrt='roll', rows=ar, cols=ar)
+        self.declare_partials(of='dx_dt', wrt='pitch', rows=ar, cols=ar)
+        self.declare_partials(of='dx_dt', wrt='yaw', rows=ar, cols=ar)
+
+        self.declare_partials(of='dy_dt', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='dy_dt', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='dy_dt', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='dy_dt', wrt='roll', rows=ar, cols=ar)
+        self.declare_partials(of='dy_dt', wrt='pitch', rows=ar, cols=ar)
+        self.declare_partials(of='dy_dt', wrt='yaw', rows=ar, cols=ar)
+
+        self.declare_partials(of='dz_dt', wrt='u', rows=ar, cols=ar)
+        self.declare_partials(of='dz_dt', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='dz_dt', wrt='w', rows=ar, cols=ar)
+        self.declare_partials(of='dz_dt', wrt='roll', rows=ar, cols=ar)
+        self.declare_partials(of='dz_dt', wrt='pitch', rows=ar, cols=ar)
+        
+
+    def compute(self, inputs, outputs):
+        """
+        Compute function for EOM. 
+        TODO: Same as above, potentially rewrite equations for \
+              matrix form, and add potential assymetry to moment \
+              of inertia matrix.
+
+        """
+
+        # inputs
+
+        mass = inputs['mass']
+        u = inputs['u'] # u
+        v = inputs['v'] # v
+        w = inputs['w'] # w
+        roll_ang_vel = inputs['roll_ang_vel'] # p
+        pitch_ang_vel = inputs['pitch_ang_vel'] # q
+        yaw_ang_vel = inputs['yaw_ang_vel'] # r
+        roll = inputs['roll'] # phi
+        pitch = inputs['pitch'] # theta
+        yaw = inputs['yaw'] # psi
+        x = inputs['x'] # p1
+        y = inputs['y'] # p2
+        z = inputs['z'] # p3
+        # time = inputs['time']
+        g = inputs['g']
+        Fx = inputs['Fx']
+        Fy = inputs['Fy']
+        Fz = inputs['Fz']
+        lx = inputs['lx'] # l
+        ly = inputs['ly'] # m
+        lz = inputs['lz'] # n
+        J_xz = inputs['J_xz']
+        J_xx = inputs['J_xx']
+        J_yy = inputs['J_yy']
+        J_zz = inputs['J_zz']
+
+        # Resolve gravity in body coordinate system -- denoted with subscript 'b'
+        gx_b = -np.sin(pitch) * g
+        gy_b = np.sin(roll) * np.cos(pitch) * g
+        gz_b = np.cos(roll) * np.cos(pitch) * g
+
+        # TODO: could add external forces and moments here if needed
+
+        # Denominator for roll and yaw rate equations
+        Den = J_xx * J_zz - J_xz**2
+
+        # roll-axis velocity equation
+
+        dx_accel = 1 / mass * Fx + gx_b - w * pitch_ang_vel + v * yaw_ang_vel
+
+        # pitch-axis velocity equation
+
+        dy_accel = 1 / mass * Fy + gy_b - u * yaw_ang_vel + w * roll_ang_vel
+
+        # yaw-axis velocity equation
+
+        dz_accel = 1 / mass * Fz + gz_b - v * roll_ang_vel + u * pitch_ang_vel
+
+        # Roll equation
+
+        roll_accel = (J_xz * (J_xx - J_yy + J_zz) * roll_ang_vel * pitch_ang_vel - 
+                   (J_zz * (J_zz - J_yy) + J_xz**2) * pitch_ang_vel * yaw_ang_vel + 
+                   J_zz * lx + 
+                   J_xz * lz) / Den
+        
+        # Pitch equation
+
+        pitch_accel = ((J_zz - J_xx) * roll_ang_vel * yaw_ang_vel - 
+                    J_xz * (roll_ang_vel**2 - yaw_ang_vel**2) + ly) / J_yy
+        
+        # Yaw equation
+
+        yaw_accel = ((J_xx * (J_xx - J_yy) + J_xz**2) * roll_ang_vel * pitch_ang_vel + 
+                  J_xz * (J_xx - J_yy + J_zz) * pitch_ang_vel * yaw_ang_vel + 
+                  J_xz * lx + 
+                  J_xz * lz) / Den
+        
+        # Kinematic equations
+        
+        roll_angle_rate_eq = roll_ang_vel + np.sin(roll) * np.tan(pitch) * pitch_ang_vel + \
+                             np.cos(roll) * np.tan(pitch) * yaw_ang_vel
+        
+        pitch_angle_rate_eq = np.cos(roll) * pitch_ang_vel - np.sin(roll) * yaw_ang_vel
+
+        yaw_angle_rate_eq = np.sin(roll) / np.cos(pitch) * pitch_ang_vel + \
+                            np.cos(roll) / np.cos(pitch) * yaw_ang_vel
+
+        # Position equations
+
+        dx_dt = np.cos(pitch) * np.cos(yaw) * u + \
+                (-np.cos(roll) * np.sin(yaw) + np.sin(roll) * np.sin(pitch) * np.cos(yaw)) * v + \
+                (np.sin(roll) * np.sin(yaw) + np.cos(roll) * np.sin(pitch) * np.cos(yaw)) * w
+        
+        dy_dt = np.cos(pitch) * np.sin(yaw) * u + \
+                (np.cos(roll) * np.cos(yaw) + np.sin(roll) * np.sin(pitch) * np.sin(yaw)) * v + \
+                (-np.sin(roll) * np.cos(yaw) + np.cos(roll) * np.sin(pitch) * np.sin(yaw)) * w
+        
+        dz_dt = -np.sin(pitch) * u + \
+                np.sin(roll) * np.cos(pitch) * v + \
+                np.cos(roll) * np.cos(pitch) * w
+
+        outputs['dx_accel'] = dx_accel
+        outputs['dy_accel'] = dy_accel
+        outputs['dz_accel'] = dz_accel
+        outputs['roll_accel'] = roll_accel
+        outputs['pitch_accel'] = pitch_accel
+        outputs['yaw_accel'] = yaw_accel
+        outputs['roll_angle_rate_eq'] = roll_angle_rate_eq
+        outputs['pitch_angle_rate_eq'] = pitch_angle_rate_eq
+        outputs['yaw_angle_rate_eq'] = yaw_angle_rate_eq
+        outputs['dx_dt'] = dx_dt
+        outputs['dy_dt'] = dy_dt
+        outputs['dz_dt'] = dz_dt
+        
+    
+    def compute_partials(self, inputs, J):
+        mass = inputs['mass']
+        u = inputs['u'] # u
+        v = inputs['v'] # v
+        w = inputs['w'] # w
+        roll_ang_vel = inputs['roll_ang_vel'] # p
+        pitch_ang_vel = inputs['pitch_ang_vel'] # q
+        yaw_ang_vel = inputs['yaw_ang_vel'] # r
+        roll = inputs['roll'] # phi
+        pitch = inputs['pitch'] # theta
+        yaw = inputs['yaw'] # psi
+        x = inputs['x'] # p1
+        y = inputs['y'] # p2
+        z = inputs['z'] # p3
+        # time = inputs['time']
+        g = inputs['g']
+        Fx = inputs['Fx']
+        Fy = inputs['Fy']
+        Fz = inputs['Fz']
+        lx = inputs['lx'] # l
+        ly = inputs['ly'] # m
+        lz = inputs['lz'] # n
+        J_xz = inputs['J_xz']
+        J_xx = inputs['J_xx']
+        J_yy = inputs['J_yy']
+        J_zz = inputs['J_zz']
+
+        # for roll and yaw
+        Den = J_xx * J_zz - J_xz**2
+
+        J['dx_accel', 'mass'] = -Fx / mass**2
+        J['dx_accel', 'Fx'] = 1 / mass
+        J['dx_accel', 'v'] = yaw_ang_vel
+        J['dx_accel', 'w'] = -pitch_ang_vel
+        J['dx_accel', 'yaw_ang_vel'] = v
+        J['dx_accel', 'pitch_ang_vel'] = -w
+        J['dx_accel', 'g'] = -np.sin(pitch)
+        J['dx_accel', 'pitch'] = -np.cos(pitch) * g
+
+        J['dy_accel', 'mass'] = -Fy / mass**2
+        J['dy_accel', 'Fy'] = 1 / mass
+        J['dy_accel', 'u'] = -yaw_ang_vel
+        J['dy_accel', 'w'] = roll_ang_vel
+        J['dy_accel', 'yaw_ang_vel'] = -u
+        J['dy_accel', 'roll_ang_vel'] = w
+        J['dy_accel', 'g'] = np.sin(roll) * np.cos(pitch)
+        J['dy_accel', 'roll'] = np.cos(roll) * np.cos(pitch) * g
+        J['dy_accel', 'pitch'] = -np.sin(roll) * np.sin(pitch) * g
+
+        J['dz_accel', 'mass'] = -Fz / mass**2
+        J['dz_accel', 'Fz'] = 1 / mass
+        J['dz_accel', 'v'] = -roll_ang_vel
+        J['dz_accel', 'u'] = pitch_ang_vel
+        J['dz_accel', 'roll_ang_vel'] = -v
+        J['dz_accel', 'pitch_ang_vel'] = u
+        J['dz_accel', 'g'] = np.cos(roll) * np.cos(pitch)
+        J['dz_accel', 'roll'] = -np.sin(roll) * np.cos(pitch) * g
+        J['dz_accel', 'pitch'] = -np.cos(roll) * np.sin(pitch) * g
+
+        J['roll_accel', 'J_xz'] = (Den * (
+            (J_xx - J_yy + J_zz) * roll_ang_vel * pitch_ang_vel - 2 * J_xz * pitch_ang_vel * yaw_ang_vel + lz) - (
+                J_xz * (J_xx - J_yy + J_zz) * roll_ang_vel * pitch_ang_vel - 
+                   (J_zz * (J_zz - J_yy) + J_xz**2) * pitch_ang_vel * yaw_ang_vel + 
+                   J_zz * lx + 
+                   J_xz * lz) * -2 * J_xz) / Den**2
+        J['roll_accel', 'J_xx'] = (Den * (
+            J_xz * roll_ang_vel * pitch_ang_vel
+        ) - (J_xz * (J_xx - J_yy + J_zz) * roll_ang_vel * pitch_ang_vel - 
+                   (J_zz * (J_zz - J_yy) + J_xz**2) * pitch_ang_vel * yaw_ang_vel + 
+                   J_zz * lx + 
+                   J_xz * lz) * J_zz) / Den**2
+        J['roll_accel', 'J_yy'] = (-J_xz * roll_ang_vel * pitch_ang_vel + J_zz * pitch_ang_vel * yaw_ang_vel) / Den
+        J['roll_accel', 'J_zz'] = (Den * (
+            J_xz * roll_ang_vel * pitch_ang_vel - 2 * J_zz * pitch_ang_vel * yaw_ang_vel + J_yy * pitch_ang_vel * yaw_ang_vel + lx
+        ) - (J_xz * (J_xx - J_yy + J_zz) * roll_ang_vel * pitch_ang_vel - 
+                   (J_zz * (J_zz - J_yy) + J_xz**2) * pitch_ang_vel * yaw_ang_vel + 
+                   J_zz * lx + 
+                   J_xz * lz) * J_xx) / Den**2
+        J['roll_accel', 'roll_ang_vel'] = (J_xz * (J_xx - J_yy + J_zz) * pitch_ang_vel) / Den
+        J['roll_accel', 'pitch_ang_vel'] = (J_xz * (J_xx - J_yy + J_zz) * roll_ang_vel - (J_zz * (J_zz - J_yy) + J_xz**2) * yaw_ang_vel) / Den
+        J['roll_accel', 'yaw_ang_vel'] = -((J_zz * (J_zz - J_yy) + J_xz**2) * pitch_ang_vel) / Den
+        J['roll_accel', 'lx'] = J_zz / Den
+        J['roll_accel', 'lz'] = J_xz / Den
+
+        J['pitch_accel', 'J_xz'] = -(roll_ang_vel**2 - yaw_ang_vel**2) / J_yy
+        J['pitch_accel', 'J_xx'] = -(roll_ang_vel * yaw_ang_vel) / J_yy
+        J['pitch_accel', 'J_yy'] = -((J_zz - J_xx) * roll_ang_vel * yaw_ang_vel - 
+                    J_xz * (roll_ang_vel**2 - yaw_ang_vel**2) + ly) / J_yy**2
+        J['pitch_accel', 'J_zz'] = roll_ang_vel * yaw_ang_vel / J_yy
+        J['pitch_accel', 'roll_ang_vel'] = ((J_zz - J_xx) * yaw_ang_vel - 2 * J_xz * roll_ang_vel) / J_yy
+        J['pitch_accel', 'yaw_ang_vel'] = ((J_zz - J_xx) * roll_ang_vel + 2 * J_xz * yaw_ang_vel) / J_yy
+        J['pitch_accel', 'ly'] = 1 / J_yy
+
+        J['yaw_accel', 'J_xz'] = (Den * (
+            2 * J_xz * roll_ang_vel * pitch_ang_vel + (J_xx - J_yy + J_zz) * pitch_ang_vel * yaw_ang_vel + lx + lz
+        ) - ((J_xx * (J_xx - J_yy) + J_xz**2) * roll_ang_vel * pitch_ang_vel + 
+                  J_xz * (J_xx - J_yy + J_zz) * pitch_ang_vel * yaw_ang_vel + 
+                  J_xz * lx + 
+                  J_xz * lz) * -2 * J_xz) / Den**2
+        J['yaw_accel', 'J_xx'] = (Den * (
+            2 * J_xx * roll_ang_vel * pitch_ang_vel - J_yy * roll_ang_vel * pitch_ang_vel + J_xz * pitch_ang_vel * yaw_ang_vel
+        ) - ((J_xx * (J_xx - J_yy) + J_xz**2) * roll_ang_vel * pitch_ang_vel + 
+                  J_xz * (J_xx - J_yy + J_zz) * pitch_ang_vel * yaw_ang_vel + 
+                  J_xz * lx + 
+                  J_xz * lz) * J_zz) / Den**2
+        J['yaw_accel', 'J_yy'] = (-J_xx * roll_ang_vel * pitch_ang_vel - J_xz * pitch_ang_vel * yaw_ang_vel) / Den
+        J['yaw_accel', 'J_zz'] = (Den * (
+            J_xz * pitch_ang_vel * yaw_ang_vel
+        ) - ((J_xx * (J_xx - J_yy) + J_xz**2) * roll_ang_vel * pitch_ang_vel + 
+                  J_xz * (J_xx - J_yy + J_zz) * pitch_ang_vel * yaw_ang_vel + 
+                  J_xz * lx + 
+                  J_xz * lz) * J_xx) / Den**2
+        J['yaw_accel', 'roll_ang_vel'] = ((J_xx * (J_xx - J_yy) + J_xz**2) * pitch_ang_vel) / Den
+        J['yaw_accel', 'pitch_ang_vel'] = ((J_xx * (J_xx - J_yy) + J_xz**2) * roll_ang_vel + J_xz * (J_xx - J_yy + J_zz) * yaw_ang_vel) / Den
+        J['yaw_accel', 'yaw_ang_vel'] = (J_xz * (J_xx - J_yy + J_zz) * pitch_ang_vel) / Den
+        J['yaw_accel', 'lx'] = J_xz / Den
+        J['yaw_accel', 'lz'] = J_xz / Den
+
+        J['roll_angle_rate_eq', 'roll_ang_vel'] = 1 
+        J['roll_angle_rate_eq', 'pitch_ang_vel'] = np.sin(roll) * np.tan(pitch)
+        J['roll_angle_rate_eq', 'yaw_ang_vel'] = np.cos(roll) * np.tan(pitch)
+        J['roll_angle_rate_eq', 'roll'] = np.cos(roll) * np.tan(pitch) * pitch_ang_vel - np.sin(roll) * np.tan(pitch) * yaw_ang_vel
+        J['roll_angle_rate_eq', 'pitch'] = np.sin(roll) * (1 / np.cos(pitch)**2) * pitch_ang_vel + np.cos(roll) * (1 / np.cos(pitch)**2) * yaw_ang_vel
+        
+        J['pitch_angle_rate_eq', 'pitch_ang_vel'] = np.cos(roll)
+        J['pitch_angle_rate_eq', 'yaw_ang_vel'] = -np.sin(roll)
+        J['pitch_angle_rate_eq', 'roll'] = -np.sin(roll) * pitch_ang_vel - np.cos(roll) * yaw_ang_vel
+
+        J['yaw_angle_rate_eq', 'pitch_ang_vel'] = np.sin(roll) / np.cos(pitch)
+        J['yaw_angle_rate_eq', 'yaw_ang_vel'] = np.cos(roll) / np.cos(pitch)
+        J['yaw_angle_rate_eq', 'roll'] = np.cos(roll) / np.cos(pitch) * pitch_ang_vel - np.sin(roll) / np.cos(pitch) * yaw_ang_vel
+        J['yaw_angle_rate_eq', 'pitch'] = np.sin(roll) * (np.tan(pitch) / np.cos(pitch)) * pitch_ang_vel + np.cos(roll) * (np.tan(pitch) / np.cos(pitch)) * yaw_ang_vel
+
+        # note: d/dx tan(x) = sec^2(x) = 1 / cos^2(x)
+        # note: d/dx 1 / cos(x) = d/dx sec(x) = sec(x)tan(x) = tan(x) / cos(x)
+
+        J['dx_dt', 'u'] = np.cos(pitch) * np.cos(yaw)
+        J['dx_dt', 'v'] = -np.cos(roll) * np.sin(yaw) + np.sin(roll) * np.sin(pitch) * np.cos(yaw)
+        J['dx_dt', 'w'] = np.sin(roll) * np.sin(yaw) + np.cos(roll) * np.sin(pitch) * np.cos(yaw)
+        J['dx_dt', 'roll'] = (np.sin(roll) * np.sin(yaw) + np.cos(roll) * np.sin(pitch) * np.cos(yaw)) * v + \
+                             (np.cos(roll) * np.sin(yaw) - np.sin(roll) * np.sin(pitch) * np.cos(yaw)) * w
+        J['dx_dt', 'pitch'] = -np.sin(pitch) * np.cos(yaw) * u + \
+                              np.sin(roll) * np.cos(pitch) * np.cos(yaw) * v + \
+                              np.cos(roll) * np.cos(pitch) * np.cos(yaw) * w
+        J['dx_dt', 'yaw'] = -np.cos(pitch) * np.sin(yaw) * u + \
+                            (-np.cos(roll) * np.cos(yaw) - np.sin(roll) * np.sin(pitch) * np.sin(yaw)) * v + \
+                            (np.sin(roll) * np.cos(yaw) - np.cos(roll) * np.sin(pitch) * np.sin(yaw)) * w
+        
+        J['dy_dt', 'u'] = np.cos(pitch) * np.sin(yaw)
+        J['dy_dt', 'v'] = np.cos(roll) * np.cos(yaw) + np.sin(roll) * np.sin(pitch) * np.sin(yaw)
+        J['dy_dt', 'w'] = -np.sin(roll) * np.cos(yaw) + np.cos(roll) * np.sin(pitch) * np.sin(yaw)
+        J['dy_dt', 'roll'] = (-np.sin(roll) * np.cos(yaw) + np.cos(roll) * np.sin(pitch) * np.sin(yaw)) * v + \
+                             (-np.cos(roll) * np.cos(yaw) - np.sin(roll) * np.sin(pitch) * np.sin(yaw)) * w
+        J['dy_dt', 'pitch'] = -np.sin(pitch) * np.sin(yaw) * u + \
+                              np.sin(roll) * np.cos(pitch) * np.sin(yaw) * v + \
+                              np.cos(roll) * np.cos(pitch) * np.sin(yaw) * w
+        J['dy_dt', 'yaw'] = np.cos(pitch) * np.cos(yaw) * u + \
+                            (-np.cos(roll) * np.sin(yaw) + np.sin(roll) * np.sin(pitch) * np.cos(yaw)) * v + \
+                            (np.sin(roll) * np.sin(yaw) + np.cos(roll) * np.sin(pitch) * np.cos(yaw)) * w
+        
+        J['dz_dt', 'u'] = -np.sin(pitch)
+        J['dz_dt', 'v'] = np.sin(roll) * np.cos(pitch)
+        J['dz_dt', 'w'] = np.cos(roll) * np.cos(pitch)
+        J['dz_dt', 'roll'] = np.cos(roll) * np.cos(pitch) * v - \
+                             np.sin(roll) * np.cos(pitch) * w
+        J['dz_dt', 'pitch'] = -np.cos(pitch) * u - \
+                              np.sin(roll) * np.sin(pitch) * v - \
+                              np.cos(roll) * np.sin(pitch) * w
+                             
+
+
+
+
+if __name__ == "__main__":
+
+    p = om.Problem()
+    p.model = om.Group()
+    des_vars = p.model.add_subsystem('des_vars', om.IndepVarComp(), promotes=['*'])
+
+    des_vars.add_output('mass', 3.0, units='kg')
+    des_vars.add_output('u', 0.1, units='m/s')
+    des_vars.add_output('v', 0.7, units='m/s')
+    des_vars.add_output('w', 0.12, units='m/s')
+    des_vars.add_output('roll_ang_vel', 0.1, units='rad/s')
+    des_vars.add_output('pitch_ang_vel', 0.9, units='rad/s')
+    des_vars.add_output('yaw_ang_vel', 0.12, units='rad/s')
+    des_vars.add_output('roll', 0.9, units='rad')
+    des_vars.add_output('pitch', 0.19, units='rad')
+    des_vars.add_output('yaw', 0.70, units='rad')
+    des_vars.add_output('g', 9.81, units='m/s**2')
+    des_vars.add_output('Fx', 0.1, units='N')
+    des_vars.add_output('Fy', 0.9, units='N')
+    des_vars.add_output('Fz', 0.12, units='N')
+    des_vars.add_output('lx', 3.0, units='N*m')
+    des_vars.add_output('ly', 4.0, units='N*m')
+    des_vars.add_output('lz', 5.0, units='N*m')
+    des_vars.add_output('J_xz', 9.0, units='kg*m**2')
+    des_vars.add_output('J_xx', 50.0, units='kg*m**2')
+    des_vars.add_output('J_yy', 51.0, units='kg*m**2')
+    des_vars.add_output('J_zz', 52.0, units='kg*m**2')
+
+    p.model.add_subsystem('SixDOF_EOM', SixDOF_EOM(num_nodes=1), promotes=['*'])
+
+    p.setup(check=False, force_alloc_complex=True)
+
+    p.run_model()
+
+    dx_accel = p.get_val('dx_accel')
+    dy_accel = p.get_val('dy_accel')
+    dz_accel = p.get_val('dz_accel')
+    roll_accel = p.get_val('roll_accel')
+    pitch_accel = p.get_val('pitch_accel')
+    yaw_accel = p.get_val('yaw_accel')
+    roll_angle_rate_eq = p.get_val('roll_angle_rate_eq')
+    pitch_angle_rate_eq = p.get_val('pitch_angle_rate_eq')
+    yaw_angle_rate_eq = p.get_val('yaw_angle_rate_eq')
+    dx_dt = p.get_val('dx_dt')
+    dy_dt = p.get_val('dy_dt')
+    dz_dt = p.get_val('dz_dt')
+
+    print(f"Accelerations in x,y,z: {dx_accel}, {dy_accel}, {dz_accel}")
+    print(f"Euler angle accels in roll, pitch, yaw: {roll_accel}, {pitch_accel}, {yaw_accel}")
+    print(f"Euler angular rates: {roll_angle_rate_eq}, {pitch_angle_rate_eq}, {yaw_angle_rate_eq}")
+    print(f"velocities: {dx_dt}, {dy_dt}, {dz_dt}")
+
+    p.check_partials(compact_print=True, show_only_incorrect=True, method='cs')
+
+
diff --git a/aviary/mission/sixdof/six_dof_ODE.py b/aviary/mission/sixdof/six_dof_ODE.py
new file mode 100644
index 000000000..a3e12f60e
--- /dev/null
+++ b/aviary/mission/sixdof/six_dof_ODE.py
@@ -0,0 +1,171 @@
+import numpy as np
+import openmdao.api as om
+
+from aviary.mission.base_ode import BaseODE as _BaseODE
+
+from aviary.variable_info.enums import AnalysisScheme, SpeedType, ThrottleAllocation
+from aviary.variable_info.variables import Aircraft, Dynamic, Mission
+from aviary.constants import GRAV_METRIC_FLOPS
+
+from aviary.mission.sixdof.six_dof_EOM import SixDOF_EOM
+from aviary.mission.sixdof.force_component_calc import ForceComponentResolver
+
+class SixDOF_ODE(_BaseODE):
+
+    def initialize(self):
+        super().initialize()
+
+    def setup(self):
+        options = self.options
+        nn = options['num_nodes']
+        analysis_scheme = options['analysis_scheme']
+
+        self.add_subsystem(
+            'true_airspeed_comp',
+            subsys=om.ExecComp(
+                'true_airspeed = (u**2 + v**2 + w**2)**0.5',
+                true_airspeed = {'units': 'm/s', 'shape': (nn,)},
+                u = {'units': 'm/s', 'shape': (nn,)},
+                v = {'units': 'm/s', 'shape': (nn,)}, 
+                w = {'units': 'm/s', 'shape': (nn,)},
+                has_diag_partials=True
+            ),
+            promotes_inputs=[
+                'u',
+                'v',
+                'w',
+            ],
+            promotes_outputs=[
+                'true_airspeed'
+            ]
+        )
+
+        self.add_atmosphere(input_speed_type=SpeedType.TAS) 
+
+        sub1 = self.add_subsystem(
+            'solver_sub',
+            om.Group(),
+            promotes=['*']
+        )
+
+        self.add_subsystem(
+            'thrust_comp',
+            subsys=om.ExecComp(
+            'T_z = cos(roll) * cos(pitch) * mass * g',
+            roll = {'units': 'rad', 'shape': (nn,)},
+            pitch = {'units': 'rad', 'shape': (nn,)},
+            mass = {'units': 'kg', 'shape': (nn,)},
+            g = {'val': GRAV_METRIC_FLOPS, 'units': 'm/s**2'},
+            promotes_inputs=[
+                ('roll', 'roll'),
+                ('pitch', 'pitch'),
+                ('mass', 'mass'),
+            ],
+            promotes_outputs=[
+                ('T_z', 'T_z'),
+            ]
+        )) # body-relative CS
+
+        comp = om.BalanceComp(
+            name=Dynamic.Vehicle.Propulsion.THROTTLE,
+            units='unitless',
+            val=np.ones((nn,)),
+            lhs_name='T_z',
+            rhs_name= Dynamic.Vehicle.Propulsion.THRUST_TOTAL,
+            eq_units='N',
+            normalize=False,
+        )
+
+        sub1.add_subsystem(
+            'throttle_balance',
+            subsys=comp,
+            promotes_inputs=['*'],
+            promotes_outputs=['*'],
+        )
+
+        self.add_core_subsystems(solver_group=sub1)
+
+        self.add_external_subsystems(solver_group=sub1)
+
+        self.add_subsystem(
+            'sum_forces_comp',
+            ForceComponentResolver(num_nodes=nn),
+            promotes_inputs=[
+                'u',
+                'v',
+                'w',
+                'drag',
+                'lift',
+                'side',
+                'thrust',
+                'heading_angle',
+                Dynamic.Mission.FLIGHT_PATH_ANGLE,
+            ],
+            promotes_outputs=[
+                'Fx',
+                'Fy',
+                'Fz',
+            ]
+        )
+
+        self.add_subsystem(
+            'SixDOF_EOM',
+            SixDOF_EOM(num_nodes=nn),
+            promotes_inputs=[
+                'mass',
+                'u',
+                'v',
+                'w',
+                'roll_ang_vel',
+                'pitch_ang_vel',
+                'yaw_ang_vel',
+                'roll',
+                'pitch',
+                'yaw',
+                'g',
+                'Fx',
+                'Fy',
+                'Fz',
+                'lx',
+                'ly',
+                'lz',
+                'J_xz',
+                'J_xx',
+                'J_yy',
+                'J_zz'
+            ],
+            promotes_outputs=[
+                'dx_accel',
+                'dy_accel',
+                'dz_accel',
+                'roll_accel',
+                'pitch_accel',
+                'yaw_accel',
+                'roll_angle_rate_eq',
+                'pitch_angle_rate_eq',
+                'yaw_angle_rate_eq',
+                'dx_dt',
+                'dy_dt',
+                'dz_dt'
+            ]
+        )
+
+        # self.set_input_defaults(Dynamic.Atmosphere.MACH, val=np.ones(nn), units='unitless')
+        # self.set_input_defaults(Dynamic.Vehicle.MASS, val=np.ones(nn), units='kg')
+        # self.set_input_defaults(Dynamic.Mission.VELOCITY, val=np.ones(nn), units='m/s')
+        # self.set_input_defaults(Dynamic.Mission.ALTITUDE, val=np.ones(nn), units='m')
+        # self.set_input_defaults(Dynamic.Mission.ALTITUDE_RATE, val=np.ones(nn), units='m/s')
+
+        print_level = 0 if analysis_scheme is AnalysisScheme.SHOOTING else 2
+
+        sub1.nonlinear_solver = om.NewtonSolver(
+            solve_subsystems=True,
+            atol=1.0e-10,
+            rtol=1.0e-10,
+        )
+        sub1.nonlinear_solver.linesearch = om.BoundsEnforceLS()
+        sub1.linear_solver = om.DirectSolver(assemble_jac=True)
+        sub1.nonlinear_solver.options['err_on_non_converge'] = True
+        sub1.nonlinear_solver.options['iprint'] = print_level
+
+        self.options['auto_order'] = True
\ No newline at end of file
diff --git a/aviary/mission/sixdof/test/test_six_dof_ODE.py b/aviary/mission/sixdof/test/test_six_dof_ODE.py
new file mode 100644
index 000000000..c0eecf290
--- /dev/null
+++ b/aviary/mission/sixdof/test/test_six_dof_ODE.py
@@ -0,0 +1,24 @@
+import unittest 
+import openmdao.api as om
+import numpy as np
+
+from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal
+
+from aviary.mission.sixdof.six_dof_ODE import SixDOF_ODE
+from aviary.variable_info.variables import Aircraft, Dynamic
+
+
+class SixDOFODETestCase(unittest.TestCase):
+    """
+    Test 6-degree of freedom ODE.
+
+    """
+
+    def setUp(self):
+        self.prob = om.Problem()
+
+        self.sys = self.prob.model = SixDOF_ODE(
+            num_nodes=1, 
+            
+        )
+
diff --git a/pyproject.toml b/pyproject.toml
index 86a77e3c5..7b3b31033 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,14 +7,13 @@ name = "aviary"
 dynamic = ["version"]
 readme = "README.md"
 license = "Apache-2.0"
-requires_python = ">=3.9"
 dependencies = [
     "dymos>=1.8.1",
     "hvplot",
     "importlib_resources",
     "matplotlib",
     "numpy<2",
-    "openmdao>=3.37.0",
+    "openmdao>=3.36.0",
     "pandas",
     "panel>=1.0.0",
     "parameterized",
@@ -22,19 +21,25 @@ dependencies = [
 ]
 
 [project.optional-dependencies]
-docs = [
-    "jupyter-book",
-    "itables"
-]
-dev = [
+all = [
+    "ambiance",
+    "itables",
+    "myst-nb",
+    "openaerostruct",
     "pre-commit",
+    "sphinx_book_theme==1.1.0",
     "testflo",
+]
+examples = [
     "ambiance",
+    "itables",
     "openaerostruct",
 ]
-all = [
-    "aviary[docs]",
-    "aviary[dev]",
+test = [
+    "myst-nb",
+    "pre-commit",
+    "sphinx_book_theme==1.1.0",
+    "testflo",
 ]
 
 [project.scripts]
@@ -60,7 +65,7 @@ quote-style = "single"
 [tool.ruff.lint]
 # isort, pydocstyle
 extend-select = ["I", "D"]
-# disabling these rules help current Aviary code pass a lint check
+# disabling these rules will help current Aviary code pass a pre-commit lint check
 extend-ignore = [
     "D100",
     "D101",