modify RV class to have thruster capability outside of atmosphere, fu…

…nctioning E2E reentry test with control

src/gncpy/dynamics/basic/reentry_vehicle.py

@@ -39,6 +39,7 @@ def __init__(
         origin_lat=-77.0563,
         origin_long=38.8719,
         CD0 = 0.25,
+        has_thrusters = True,
         **kwargs,
     ):
         """Initialize an object.
@@ -58,6 +59,10 @@ def __init__(
         CD0: float, optional
             zero-lift drag coefficient; only needed for determining induced drag for manuevering RV's. 
             If control input is None, this parameter affects nothing and is inconsequential. 
+        has_thrusters: bool, optional
+            Determines if the vehicle will respond to commanded turn/climb accelerations beyond what it can do with aero
+            forces alone (typically applies in upper/exo-atmosphere). If false, control commands will be limited to 
+            accelerations possible from aero forces alone. 
         **kwargs : dict
             Additional arguments for the parent class.
         """
@@ -66,9 +71,10 @@ def __init__(
         self.origin_lat = origin_lat
         self.origin_long = origin_long
         self.CD0 = CD0
+        self.has_thrusters = has_thrusters
 
         '''
-        control model is accelerations (in m/s) in the velocity-turn-climb frame:
+        control model is specific forces (in m/s2) in the velocity-turn-climb frame:
         u = np.array([a_thrust, # in direction of vehicle velocity
                       a_turn,   # in direction (left turn positive) perpendicular to velocity in the horizontal ENU plane
                       a_climb   # in direction (up positive) perpendicular to velocity and turn
@@ -107,12 +113,18 @@ def ENU_control_acceleration(t, x, u, *args):
 
             # limit maximum lift, assume CLMax = 3
             # (set high to poorly account for fact that thrust can be unaligned with velocity to help turn)
+            # If aero lift not enough to satisfy control, assume remaining maneuvering 
+            # accelerations provided by thrusters
             u_limited = copy.deepcopy(u).astype(float)
             total_lift_acceleration = np.linalg.norm(u[1:])
             lift_accel_max = 3*q*self.drag_parameter/self.CD0
+
             if total_lift_acceleration > lift_accel_max:
-                u_limited[1:] = u_limited[1:]*lift_accel_max/total_lift_acceleration
-                total_lift_acceleration = np.linalg.norm(u_limited[1:])
+                if self.has_thrusters:
+                    total_lift_acceleration = lift_accel_max
+                else:
+                    u_limited[1:] = u_limited[1:]*lift_accel_max/total_lift_acceleration
+                    total_lift_acceleration = np.linalg.norm(u_limited[1:])
 
             # calculate additional induced drag
             # following equation derived from CDL=CL**2/4 result from slender body potential flow theory (cite Cronvich Missile Aerodynamics)
@@ -124,7 +136,6 @@ def ENU_control_acceleration(t, x, u, *args):
 
             return a_control_ENU
 
-        # FIXME change these to calculate induced drag from maneuvering
         def g3(t, x, u, *args):
             if u is None:
                 return 0

test/unit/test_reentry_vehicle.py

@@ -22,7 +22,7 @@ def test_rv_prop():
                                  reentry_speed*np.sin(reentry_heading)*np.cos(reentry_angle),
                                  -reentry_speed*np.sin(reentry_angle)
                                  ])
-    state[0] = np.concatenate(([100000, 100000, 100000], reentry_velocity))
+    state[0] = np.concatenate(([300000, 300000, 100000], reentry_velocity))
 
     for kk, tt in enumerate(time[:-1]):
         state[kk + 1, :] = dynObj.propagate_state(tt, state[kk].reshape((-1, 1))).flatten()
@@ -65,7 +65,7 @@ def test_rv_prop():
         fig.axes[1].grid(True)
 
         fig.axes[2].plot(time, state[:, 5])
-        fig.axes[2].set_ylabel("N-vel (m/s)")
+        fig.axes[2].set_ylabel("U-vel (m/s)")
         fig.axes[2].set_xlabel("time (s)")
         fig.axes[2].grid(True)
 
@@ -87,7 +87,7 @@ def test_rv_prop():
         fig.axes[1].grid(True)
         fig.axes[1].set_xlabel("East position (m)")
 
-        fig.suptitle("Reentry vehicle simulation with RH 10m/s2 turn")
+        fig.suptitle("Reentry vehicle simulation without control input")
 
 
 def test_rv_control():
@@ -110,7 +110,7 @@ def test_rv_control():
                                  reentry_speed*np.sin(reentry_heading)*np.cos(reentry_angle),
                                  -reentry_speed*np.sin(reentry_angle)
                                  ])
-    state[0] = np.concatenate(([100000, 100000, 100000], reentry_velocity))
+    state[0] = np.concatenate(([300000, 300000, 100000], reentry_velocity))
 
     for kk, tt in enumerate(time[:-1]):
         state[kk + 1, :] = dynObj.propagate_state(
@@ -194,7 +194,7 @@ def test_rv_control():
         plt.close("all")
 
     test_rv_prop()
-    #test_rv_control()
+    test_rv_control()
 
     if DEBUG:
         plt.show()