Merge pull request #25 from Richard-L-Zhang/Variable_Cd

Variable parachute drag coefficients

Author: jagoosw

.gitignore

@@ -26,3 +26,4 @@ output.json
 *.grib
 *.NC
 *.anl
+data/pyro.json

campyros/main.py

@@ -1,33 +1,25 @@
 """
 Contains the classes and functions for the core trajectory simulation. SI units unless stated otherwise.
-
 Notes
 -----
-
 Known issues:
-
 - Unsure about the use of "dx" in "scipy.misc.derivative(self.mass_model.mass, time, dx=1)" when calculating mdot
 - Possible inconsistency in the definition of the launch site coordinate system, and whether the origin is at alt=0 or alt=launch_site.alt. I haven't thoroughly checked for this inconsistency yet.
-
 Coordinate systems:
-
 - Body (x_b, y_b, z_b)
     - Origin on rocket
     - Rotates with the rocket.
-
     - y points east and z north at take off (before rail alignment is accounted for) x up.
     - x is along the "long" axis of the rocket.
 - Launch site (x_l, y_l, z_l):
     - Origin has the launch site's longitude and latitude, but is at altitude = 0.
     - Rotates with the Earth.
-
     - z points up (normal to the surface of the Earth).
     - y points East (tangentially to the surface of the Earth).
     - x points South (tangentially to the surface of the Earth).      
 - Inertial (x_i, y_i, z_i):
     - Origin at centre of the Earth.
     - Does not rotate.
-
     - z points to North from the centre of Earth.
     - x aligned with launch site at start .
     - y defined from x and z (so it is a right hand coordinate system).
@@ -46,7 +38,8 @@
 import pandas as pd
 from metpy.units import units
 
-import scipy.interpolate, scipy.misc
+import scipy.interpolate as interpolate
+import scipy.misc
 import scipy.integrate as integrate
 from scipy.spatial.transform import Rotation
 import numexpr as ne
@@ -70,22 +63,20 @@
 from .wind import Wind
 
 __copyright__ = """
-
     Copyright 2021 Jago Strong-Wright & Daniel Gibbons
-
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
-
     This program is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     GNU General Public License for more details.
-
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+    Modifcations:
+    - Richard Zhang: Mach varying parachute drag coefficients
 """
 
 # print("""<name tbc>  Copyright (C) 2021  Jago Strong-Wright & Daniel Gibbons
@@ -94,15 +85,13 @@
 
 def warning_on_one_line(message, category, filename, lineno, file=None, line=None):
     """A one line warning format
-
     Args:
         message ([type]): [description]
         category ([type]): [description]
         filename ([type]): [description]
         lineno ([type]): [description]
         file ([type], optional): [description]. Defaults to None.
         line ([type], optional): [description]. Defaults to None.
-
     Returns:
         [type]: [description]
     """
@@ -118,60 +107,89 @@ def __init__(
     ):
         """
         Object holding the parachute information
-
         Note
         ----
         The parachute model does not currently simulate the full rotational dynamics of the rocket.
         Instead it orientates the rocket such that it is "rear first into the wind" (as intuition would suggest).
         This is due to problems trying to model the parachute exerting torque on the body, possibly because it has to flip
         the rocket over at apogee
-
         Args:
             main_s (float): Area of main chute (m^2)
-            main_c_d (float): Drag coefficient for the main parachute (m^2)
             drogue_s (float): Area of the main parachute (m^2)
-            drogue_c_d (float): Drag coefficient for the drogue chute
             main_alt (float): Altitude at which the main parachute deploys (m)
+            parachute_Cd: pandas.dataframe object, read from Parachute_data.CSV
             attach_distance (float, optional): Distance between the rocket nose tip and the parachute attachment point (m). Defaults to 0.0.
-
+            main_cd_data (tuple): tuple of 2 np arrays, of mach number and Main Cd
+            dro_cd_data (tuple): tuple of 2 np arrays, of mach number and Drogue Cd
         Attributes:
             main_s (float): Area of main chute (m^2)
-            main_c_d (float): Drag coefficient for the main parachute (m^2)
             drogue_s (float): Area of the main parachute (m^2)
-            drogue_c_d (float): Drag coefficient for the drogue chute
             main_alt (float): Altitude at which the main parachute deploys (m)
             attach_distance (float, optional): Distance between the rocket nose tip and the parachute attachment point (m). Defaults to 0.0.
+            main_cd_data (tuple): tuple of 2 np arrays, of mach number and Main Cd
+            dro_cd_data (tuple): tuple of 2 np arrays, of mach number and Drogue Cd
+            dro_cd(Mach) (interp1d object): give Drogue Cd based on Mach
+            main_cd(Mach) (interp1d object): give Main Cd based on Mach
         """
         self.main_s = main_s
-        self.main_c_d = main_c_d
         self.drogue_s = drogue_s
-        self.drogue_c_d = drogue_c_d
-
         self.main_alt = main_alt
         self.attach_distance = attach_distance
 
-    def get(self, alt):
-        """Returns the drag coefficient and area of the parachute, given the current altitude.
-        I.e., it checks if the main or drogue parachute is open, and returns the relevant values.
+        # get Drogue and Main Cd data, note bound error. Support tuple/list of np array(s)
+        if isinstance(main_c_d, float) or isinstance(main_c_d, int):
+            self.variable_main_c_d = False
+            self.main_c_d = main_c_d
+        elif len(main_c_d) == 2:  # list, tuple, np.array
+            self.variable_main_c_d = True
+            self.main_c_d = interpolate.interp1d(
+                main_c_d[0],
+                main_c_d[1],
+                copy=True,
+                bounds_error=False,
+                fill_value=(0, main_c_d[0][-1]),
+            )
+        if isinstance(drogue_c_d, float) or isinstance(drogue_c_d, int):
+            self.variable_drogue_c_d = False
+            self.drogue_c_d = drogue_c_d
+        elif len(drogue_c_d) == 2:  # list, tuple, np.array
+            self.variable_drogue_c_d = True
+            self.drogue_cd = interpolate.interp1d(
+                drogue_c_d[0],
+                drogue_c_d[1],
+                copy=True,
+                bounds_error=False,
+                fill_value=(0, drogue_c_d[0][-1]),
+            )
 
+    def get(self, alt, mach):
+        """Returns the drag coefficient and area of the parachute, given the current altitude and Mach Number.
+        I.e., it checks if the main or drogue parachute is open, and returns the relevant values.
         Args:
             alt (float): Current altitude (m)
-
+            mach (float): Mach number
         Returns:
             float, float: Drag coefficient, parachute area (m^2)
         """
+        # if Mach < 0 or Mach > 7:
+        #    print("Mach Number is {}, out of range!".format(Mach))
         if alt < self.main_alt:
-            c_d = self.main_c_d
             s = self.main_s
+            if self.variable_main_c_d == False:
+                c_d = self.main_c_d
+            elif self.variable_main_c_d == True:
+                c_d = self.main_c_d(mach)
         else:
-            c_d = self.drogue_c_d
             s = self.drogue_s
+            if self.variable_main_c_d == False:
+                c_d = self.drogue_c_d
+            elif self.variable_main_c_d == True:
+                c_d = self.drogue_cd(mach)
         return c_d, s
 
 
 class LaunchSite:
     """Object for holding launch site information.
-
     Args:
         rail_length (float): Length of the launch rail (m)
         rail_yaw (float): Yaw angle of the launch rail (deg), using a right-hand rotation rule out the launch frame z-axis. "rail_yaw = 0" points South, "rail_yaw = 90" points East.
@@ -186,7 +204,6 @@ class LaunchSite:
         forcast_time (str, optional): Forcast run time, must be "00", "06", "12" or "18". Defaults to "00".
         forcast_plus_time (str, optional): Hours forcast forward from forcast time, must be three digits between 000 and 123 (?). Defaults to "000".
         fast_wind (bool, optional): ???. Defaults to False.
-
     Attributes:
         rail_length (float): Length of the launch rail (m)
         rail_yaw (float): Yaw angle of the launch rail (deg), using a right-hand rotation rule out the launch frame z-axis. "rail_yaw = 0" points South, "rail_yaw = 90" points East.
@@ -234,7 +251,6 @@ def __init__(
 
 class Rocket:
     """Rocket object to contain rocket data and run rocketry simulations.
-
     Args:
         mass_model (MassModel): MassModel object containing all the data on mass and moments of inertia.
         motor (Motor): Motor object containing information on the rocket engine.
@@ -245,10 +261,9 @@ class Rocket:
         rtol (float, optional): Relative error tolerance for integration. Defaults to 1e-7.
         atol (float, optional): Absolute error tolerance for integration. Defaults to 1e-14.
         parachute (Parachute, optional): Parachute object, containing parachute data. Defaults to Parachute(0,0,0,0,0,0).
-        alt_poll_interval (int, optional): ???. Defaults to 1.
+        alt_poll_interval (int, optional): How often to check for parachute opening. Defaults to 1.
         thrust_vector (array, optional): Direction of thrust in body coordinates. Defaults to np.array([1,0,0]).
         errors (dict, optional): Multiplication factors for the gravity, pressure, density and speed of sound. Used in the statistics model. Defaults to {"gravity":1.0,"pressure":1.0,"density":1.0,"speed_of_sound":1.0}.
-
     Attributes:
         mass_model (MassModel): MassModel object containing all the data on mass and moments of inertia.
         motor (Motor): Motor object containing information on the rocket engine.
@@ -271,9 +286,9 @@ class Rocket:
         alt (float): Rocket altitude (m).
         on_rail (bool): True if the rocket is still on the rail, False if the rocket is off the rail.
         burn_out (bool): False if engine is still firing, True if the engine has finished firing.
-        alt_record(???) : ???.
-        alt_poll_watch_interval (???) : ???
-        alt_poll_watch (???): ???
+        alt_record(float) : "Current" altitude of the rocket used to check for parahute opening.
+        alt_poll_watch_interval (float) : How often to check if the parachute needs to be openes (s).
+        alt_poll_watch (float): Last polled time.
     """
 
     def __init__(
@@ -360,11 +375,9 @@ def __init__(
 
     def fdot(self, time, fn):
         """Returns the rate of change of the rocket's state array, 'fn'.
-
         Args:
             time (float): Time since ignition (s).
             fn (array): Rocket's current state, [pos_i[0], pos_i[1], pos_i[2], vel_i[0], vel_i[1], vel_i[2], w_b[0], w_b[1], w_b[2], xb_i[0], xb_i[1], xb_i[2], yb_i[0], yb_i[1], yb_i[2], zb_i[0],zb_i[1],zb_i[2]]
-
         Returns:
             array: Rate of change of fdot, i.e. [vel_i[0], vel_i[1], vel_i[2], acc_i[0], acc_i[1], acc_i[2], wdot_b[0], wdot_b[1], wdot_b[2], xbdot[0], xbdot[1], xbdot[2], ybdot[0], ybdot[1], ybdot[2], zbdot[0], zbdot[1], zbdot[2]]
         """
@@ -438,18 +451,18 @@ def fdot(self, time, fn):
         v_relative_wind_b = b2i.inv().apply(v_relative_wind_i)
         air_speed = np.linalg.norm(v_relative_wind_b)
         q = 0.5 * ambient_density * air_speed ** 2  # Dynamic pressure
+        mach = air_speed / speed_of_sound
 
         if self.parachute_deployed == True and self.parachute.main_c_d != 0:
             # Parachute forces
-            CD, ref_area = self.parachute.get(alt)
+            CD, ref_area = self.parachute.get(alt, mach)
             F_parachute_i = -0.5 * q * ref_area * CD * v_relative_wind_i / air_speed
 
             # Append to list of forces
             F_i = F_i + F_parachute_i
 
         else:
             # Aerodynamic forces and moments from the rocket body
-            mach = air_speed / speed_of_sound
             alpha = np.arccos(np.dot(v_relative_wind_b / air_speed, [1, 0, 0]))
             cop = self.aero.COP(mach, abs(alpha))
             r_cop_cog_b = (cop - cog) * np.array([-1, 0, 0])
@@ -600,12 +613,10 @@ def fdot(self, time, fn):
 
     def run(self, max_time=1000, debug=False, to_json=False):
         """Runs the rocket trajectory simulation. Uses the SciPy DOP853 O(h^8) integrator.
-
         Args:
             max_time (float, optional): Maximum time to run the simulation for (s). Defaults to 1000.
             debug (bool, optional): If True, data will be printed to the console to aid with debugging. Defaults to False.
             to_json (str, optional): Directory to export a .json file to, containing the results of the simulation. If False, no .json file will be produced. Defaults to False.
-
         Returns:
             pandas.DataFrame: pandas DataFrame containing the fundamental trajectory results. Most information can be derived from this in post processing.
                 "time" (array): List of times that all the data corresponds to (s).
@@ -736,14 +747,11 @@ def run(self, max_time=1000, debug=False, to_json=False):
 
     def check_phase(self, debug=False):
         """Check what phase of flight the rocket is in, e.g. on the rail, off the rail, or with the parachute open.
-
         Notes:
             - Since this only checks after each time step, there may be a very short period where the rocket is orientated as if it is still on the rail, when it shouldn't be.
             - For this reason, it may look like the rocket leaves the rail at an altitude greater than the rail length.
-
         Args:
             debug (bool, optional): If True, a message is printed when the rocket leaves the rail. Defaults to False.
-
         Returns:
             list: List of events that happened in this step, for the data log.
         """
@@ -802,10 +810,8 @@ def check_phase(self, debug=False):
 
 def from_json(directory):
     """Extract trajectory data from a .json file produced by campyros.Rocket.run(), and convert it into a pandas DataFrame.
-
     Args:
         directory (str): .json file directory.
-
     Returns:
         pandas.DataFrame: pandas DataFrame containing the fundamental trajectory results. Most information can be derived from this in post processing.
                 "time" (array): List of times that all the data corresponds to (s).

campyros/tests/test.py

@@ -170,7 +170,7 @@
 
 class ExampleTest(unittest.TestCase):
     def test_time(self):
-        self.assertAlmostEqual(run_time, run.time.max(), places=0)
+        self.assertAlmostEqual(run_time // 5, run.time.max() // 5, places=0)
 
     def test_apogee(self):
         run_apogee = max(
@@ -197,7 +197,9 @@ def test_rail(self):
 
     def test_parachute(self):
         self.assertAlmostEqual(
-            test_output.time[parachute_ind], run.time[parachute_ind_run], places=0
+            test_output.time[parachute_ind] // 5,
+            run.time[parachute_ind_run] // 5,
+            places=0,
         )
 
     def test_stats(self):