diff --git a/.coveragerc b/.coveragerc new file mode 100644 index 000000000..e4405ddbc --- /dev/null +++ b/.coveragerc @@ -0,0 +1,4 @@ +[report] +exclude_also= + ; Don't complain about exceptions or warnings not being covered by tests + warnings.warn* \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index e722d335f..6c0887acd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -38,8 +38,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ### Changed - MNT: Modularize Rocket Draw [#580](https://github.com/RocketPy-Team/RocketPy/pull/580) +- DOC: Improvements of Environment docstring phrasing [#565](https://github.com/RocketPy-Team/RocketPy/pull/565) - MNT: Refactor flight prints module [#579](https://github.com/RocketPy-Team/RocketPy/pull/579) - DOC: Convert CompareFlights example notebooks to .rst files [#576](https://github.com/RocketPy-Team/RocketPy/pull/576) +- MNT: Refactor inertia calculations using parallel axis theorem [#573] (https://github.com/RocketPy-Team/RocketPy/pull/573) - ENH: Optional argument to show the plot in Function.compare_plots [#563](https://github.com/RocketPy-Team/RocketPy/pull/563) ### Fixed diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 546d446e8..b72b2cb38 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -57,8 +57,7 @@ class Environment: Environment.datum : string The desired reference ellipsoid model, the following options are available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "SIRGAS2000", then this model will be used if the user make some - typing mistake + is "SIRGAS2000". Environment.initial_east : float Launch site East UTM coordinate Environment.initial_north : float @@ -74,7 +73,7 @@ class Environment: Launch site E/W hemisphere Environment.elevation : float Launch site elevation. - Environment.date : datetime + Environment.datetime_date : datetime Date time of launch in UTC. Environment.local_date : datetime Date time of launch in the local time zone, defined by @@ -276,49 +275,70 @@ def __init__( timezone="UTC", max_expected_height=80000.0, ): - """Initialize Environment class, saving launch rail length, - launch date, location coordinates and elevation. Note that - by default the standard atmosphere is loaded until another + """Initializes the Environment class, capturing essential parameters of + the launch site, including the launch date, geographical coordinates, + and elevation. This class is designed to calculate crucial variables + for the Flight simulation, such as atmospheric air pressure, density, + and gravitational acceleration. + + Note that the default atmospheric model is the International Standard + Atmosphere as defined by ISO 2533 unless specified otherwise in + :meth:`Environment.set_atmospheric_model`. Parameters ---------- gravity : int, float, callable, string, array, optional Surface gravitational acceleration. Positive values point the - acceleration down. If None, the Somigliana formula is used to - date : array, optional - Array of length 4, stating (year, month, day, hour (UTC)) - of rocket launch. Must be given if a Forecast, Reanalysis + acceleration down. If None, the Somigliana formula is used. + See :meth:`Environment.set_gravity_model` for more information. + date : list or tuple, optional + List or tuple of length 4, stating (year, month, day, hour) in the + time zone of the parameter ``timezone``. + Alternatively, can be a ``datetime`` object specifying launch + date and time. The dates are stored as follows: + + - :attr:`Environment.local_date`: Local time of launch in + the time zone specified by the parameter ``timezone``. + + - :attr:`Environment.datetime_date`: UTC time of launch. + + Must be given if a Forecast, Reanalysis or Ensemble, will be set as an atmospheric model. + Default is None. + See :meth:`Environment.set_date` for more information. latitude : float, optional Latitude in degrees (ranging from -90 to 90) of rocket launch location. Must be given if a Forecast, Reanalysis or Ensemble will be used as an atmospheric model or if - Open-Elevation will be used to compute elevation. + Open-Elevation will be used to compute elevation. Positive + values correspond to the North. Default value is 0, which + corresponds to the equator. longitude : float, optional - Longitude in degrees (ranging from -180 to 360) of rocket + Longitude in degrees (ranging from -180 to 180) of rocket launch location. Must be given if a Forecast, Reanalysis or Ensemble will be used as an atmospheric model or if - Open-Elevation will be used to compute elevation. + Open-Elevation will be used to compute elevation. Positive + values correspond to the East. Default value is 0, which + corresponds to the Greenwich Meridian. elevation : float, optional Elevation of launch site measured as height above sea level in meters. Alternatively, can be set as 'Open-Elevation' which uses the Open-Elevation API to find elevation data. For this option, latitude and longitude must also be specified. Default value is 0. - datum : string + datum : string, optional The desired reference ellipsoidal model, the following options are available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default - is "SIRGAS2000", then this model will be used if the user make some - typing mistake. + is "SIRGAS2000". timezone : string, optional Name of the time zone. To see all time zones, import pytz and run - print(pytz.all_timezones). Default time zone is "UTC". + ``print(pytz.all_timezones)``. Default time zone is "UTC". max_expected_height : float, optional Maximum altitude in meters to keep weather data. The altitude must be above sea level (ASL). Especially useful for visualization. Can be altered as desired by doing `max_expected_height = number`. Depending on the atmospheric model, this value may be automatically - mofified. + modified. Returns ------- @@ -396,15 +416,57 @@ def set_date(self, date, timezone="UTC"): Parameters ---------- - date : Datetime - Datetime object specifying launch date and time. + date : list, tuple, datetime + List or tuple of length 4, stating (year, month, day, hour) in the + time zone of the parameter ``timezone``. See Notes for more + information. + Alternatively, can be a ``datetime`` object specifying launch + date and time. timezone : string, optional Name of the time zone. To see all time zones, import pytz and run - print(pytz.all_timezones). Default time zone is "UTC". + ``print(pytz.all_timezones)``. Default time zone is "UTC". Returns ------- None + + Notes + ----- + - If the ``date`` is given as a list or tuple, it should be in the same + time zone as specified by the ``timezone`` parameter. This local + time will be available in the attribute :attr:`Environment.local_date` + while the UTC time will be available in the attribute + :attr:`Environment.datetime_date`. + + - If the ``date`` is given as a ``datetime`` object without a time zone, + it will be assumed to be in the same time zone as specified by the + ``timezone`` parameter. However, if the ``datetime`` object has a time + zone specified in its ``tzinfo`` attribute, the ``timezone`` + parameter will be ignored. + + Examples + -------- + + Let's set the launch date as an list: + + >>> date = [2000, 1, 1, 13] # January 1st, 2000 at 13:00 UTC+1 + >>> env = Environment() + >>> env.set_date(date, timezone="Europe/Rome") + >>> print(env.datetime_date) # Get UTC time + 2000-01-01 12:00:00+00:00 + >>> print(env.local_date) + 2000-01-01 13:00:00+01:00 + + Now let's set the launch date as a ``datetime`` object: + + >>> from datetime import datetime + >>> date = datetime(2000, 1, 1, 13, 0, 0) + >>> env = Environment() + >>> env.set_date(date, timezone="Europe/Rome") + >>> print(env.datetime_date) # Get UTC time + 2000-01-01 12:00:00+00:00 + >>> print(env.local_date) + 2000-01-01 13:00:00+01:00 """ # Store date and configure time zone self.timezone = timezone @@ -458,23 +520,66 @@ def set_location(self, latitude, longitude): self.atmospheric_model_file, self.atmospheric_model_dict ) - # Return None - - def set_gravity_model(self, gravity): - """Sets the gravity model to be used in the simulation based on the - given user input to the gravity parameter. + def set_gravity_model(self, gravity=None): + """Defines the gravity model based on the given user input to the + gravity parameter. The gravity model is responsible for computing the + gravity acceleration at a given height above sea level in meters. Parameters ---------- - gravity : None or Function source - If None, the Somigliana formula is used to compute the gravity - acceleration. Otherwise, the user can provide a Function object - representing the gravity model. + gravity : int, float, callable, string, list, optional + The gravitational acceleration in m/s² to be used in the + simulation, this value is positive when pointing downwards. + The input type can be one of the following: + + - ``int`` or ``float``: The gravity acceleration is set as a\ + constant function with respect to height; + + - ``callable``: This callable should receive the height above\ + sea level in meters and return the gravity acceleration; + + - ``list``: The datapoints should be structured as\ + ``[(h_i,g_i), ...]`` where ``h_i`` is the height above sea\ + level in meters and ``g_i`` is the gravity acceleration in m/s²; + + - ``string``: The string should correspond to a path to a CSV file\ + containing the gravity acceleration data; + + - ``None``: The Somigliana formula is used to compute the gravity\ + acceleration. + + This parameter is used as a :class:`Function` object source, check\ + out the available input types for a more detailed explanation. Returns ------- Function Function object representing the gravity model. + + Notes + ----- + This method **does not** set the gravity acceleration, it only returns + a :class:`Function` object representing the gravity model. + + Examples + -------- + Let's prepare a `Environment` object with a constant gravity + acceleration: + + >>> g_0 = 9.80665 + >>> env_cte_g = Environment(gravity=g_0) + >>> env_cte_g.gravity([0, 100, 1000]) + [9.80665, 9.80665, 9.80665] + + It's also possible to variate the gravity acceleration by defining + its function of height: + + >>> R_t = 6371000 + >>> g_func = lambda h : g_0 * (R_t / (R_t + h))**2 + >>> env_var_g = Environment(gravity=g_func) + >>> g = env_var_g.gravity(1000) + >>> print(f"{g:.6f}") + 9.803572 """ if gravity is None: return self.somigliana_gravity.set_discrete( @@ -500,7 +605,7 @@ def max_expected_height(self, value): @funcify_method("height (m)", "gravity (m/s²)") def somigliana_gravity(self, height): - """Computes the gravity acceleration with the Somigliana formula. + """Computes the gravity acceleration with the Somigliana formula [1]_. An height correction is applied to the normal gravity that is accurate for heights used in aviation. The formula is based on the WGS84 ellipsoid, but is accurate for other reference ellipsoids. @@ -514,6 +619,10 @@ def somigliana_gravity(self, height): ------- Function Function object representing the gravity model. + + References + ---------- + .. [1] https://en.wikipedia.org/wiki/Theoretical_gravity#Somigliana_equation """ a = 6378137.0 # semi_major_axis f = 1 / 298.257223563 # flattening_factor diff --git a/rocketpy/environment/environment_analysis.py b/rocketpy/environment/environment_analysis.py index c15b32551..0ed638091 100644 --- a/rocketpy/environment/environment_analysis.py +++ b/rocketpy/environment/environment_analysis.py @@ -441,7 +441,7 @@ def __localize_input_dates(self): def __find_preferred_timezone(self): if self.preferred_timezone is None: - # Use local timezone based on lat lon pair + # Use local time zone based on lat lon pair try: timezonefinder = import_optional_dependency("timezonefinder") tf = timezonefinder.TimezoneFinder() diff --git a/rocketpy/motors/hybrid_motor.py b/rocketpy/motors/hybrid_motor.py index 4b28a96d2..6f0849cd0 100644 --- a/rocketpy/motors/hybrid_motor.py +++ b/rocketpy/motors/hybrid_motor.py @@ -1,3 +1,5 @@ +from rocketpy.tools import parallel_axis_theorem_from_com + from ..mathutils.function import Function, funcify_method, reset_funcified_methods from ..plots.hybrid_motor_plots import _HybridMotorPlots from ..prints.hybrid_motor_prints import _HybridMotorPrints @@ -455,22 +457,21 @@ def propellant_I_11(self): ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - solid_correction = ( - self.solid.propellant_mass - * (self.solid.center_of_propellant_mass - self.center_of_propellant_mass) - ** 2 - ) - liquid_correction = ( - self.liquid.propellant_mass - * (self.liquid.center_of_propellant_mass - self.center_of_propellant_mass) - ** 2 - ) - I_11 = ( - self.solid.propellant_I_11 - + solid_correction - + self.liquid.propellant_I_11 - + liquid_correction + solid_mass = self.solid.propellant_mass + liquid_mass = self.liquid.propellant_mass + + cm = self.center_of_propellant_mass + solid_cm_to_cm = self.solid.center_of_propellant_mass - cm + liquid_cm_to_cm = self.liquid.center_of_propellant_mass - cm + + solid_prop_inertia = self.solid.propellant_I_11 + liquid_prop_inertia = self.liquid.propellant_I_11 + + I_11 = parallel_axis_theorem_from_com( + solid_prop_inertia, solid_mass, solid_cm_to_cm + ) + parallel_axis_theorem_from_com( + liquid_prop_inertia, liquid_mass, liquid_cm_to_cm ) return I_11 diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py index ac525b379..01f728473 100644 --- a/rocketpy/motors/liquid_motor.py +++ b/rocketpy/motors/liquid_motor.py @@ -7,6 +7,7 @@ funcify_method, reset_funcified_methods, ) +from rocketpy.tools import parallel_axis_theorem_from_com from ..plots.liquid_motor_plots import _LiquidMotorPlots from ..prints.liquid_motor_prints import _LiquidMotorPrints @@ -388,10 +389,9 @@ def propellant_I_11(self): for positioned_tank in self.positioned_tanks: tank = positioned_tank.get("tank") tank_position = positioned_tank.get("position") - I_11 += ( - tank.inertia - + tank.fluid_mass - * (tank_position + tank.center_of_mass - center_of_mass) ** 2 + distance = tank_position + tank.center_of_mass - center_of_mass + I_11 += parallel_axis_theorem_from_com( + tank.inertia, tank.fluid_mass, distance ) return I_11 diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py index 6c2242a9f..3834f4a15 100644 --- a/rocketpy/motors/motor.py +++ b/rocketpy/motors/motor.py @@ -7,7 +7,7 @@ from ..mathutils.function import Function, funcify_method from ..plots.motor_plots import _MotorPlots from ..prints.motor_prints import _MotorPrints -from ..tools import tuple_handler +from ..tools import parallel_axis_theorem_from_com, tuple_handler try: from functools import cached_property @@ -513,25 +513,19 @@ def I_11(self): ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ - # Propellant inertia tensor 11 component wrt propellant center of mass - propellant_I_11 = self.propellant_I_11 - # Dry inertia tensor 11 component wrt dry center of mass + prop_I_11 = self.propellant_I_11 dry_I_11 = self.dry_I_11 - # Steiner theorem the get inertia wrt motor center of mass - propellant_I_11 += ( - self.propellant_mass - * (self.center_of_propellant_mass - self.center_of_mass) ** 2 - ) + prop_to_cm = self.center_of_propellant_mass - self.center_of_mass + dry_to_cm = self.center_of_dry_mass_position - self.center_of_mass - dry_I_11 += ( - self.dry_mass - * (self.center_of_dry_mass_position - self.center_of_mass) ** 2 + prop_I_11 = parallel_axis_theorem_from_com( + prop_I_11, self.propellant_mass, prop_to_cm ) + dry_I_11 = parallel_axis_theorem_from_com(dry_I_11, self.dry_mass, dry_to_cm) - # Sum of inertia components - return propellant_I_11 + dry_I_11 + return prop_I_11 + dry_I_11 @funcify_method("Time (s)", "Inertia I_22 (kg m²)") def I_22(self): diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 8154f2967..aaf585091 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -18,6 +18,7 @@ ) from rocketpy.rocket.components import Components from rocketpy.rocket.parachute import Parachute +from rocketpy.tools import parallel_axis_theorem_from_com class Rocket: @@ -619,6 +620,10 @@ def evaluate_dry_inertias(self): ---------- .. [1] https://en.wikipedia.org/wiki/Moment_of_inertia#Inertia_tensor """ + # Get masses + motor_dry_mass = self.motor.dry_mass + mass = self.mass + # Compute axes distances noMCM_to_CDM = ( self.center_of_mass_without_motor - self.center_of_dry_mass_position @@ -628,18 +633,18 @@ def evaluate_dry_inertias(self): ) # Compute dry inertias - self.dry_I_11 = ( - self.I_11_without_motor - + self.mass * noMCM_to_CDM**2 - + self.motor.dry_I_11 - + self.motor.dry_mass * motorCDM_to_CDM**2 + self.dry_I_11 = parallel_axis_theorem_from_com( + self.I_11_without_motor, mass, noMCM_to_CDM + ) + parallel_axis_theorem_from_com( + self.motor.dry_I_11, motor_dry_mass, motorCDM_to_CDM ) - self.dry_I_22 = ( - self.I_22_without_motor - + self.mass * noMCM_to_CDM**2 - + self.motor.dry_I_22 - + self.motor.dry_mass * motorCDM_to_CDM**2 + + self.dry_I_22 = parallel_axis_theorem_from_com( + self.I_22_without_motor, mass, noMCM_to_CDM + ) + parallel_axis_theorem_from_com( + self.motor.dry_I_22, motor_dry_mass, motorCDM_to_CDM ) + self.dry_I_33 = self.I_33_without_motor + self.motor.dry_I_33 self.dry_I_12 = self.I_12_without_motor + self.motor.dry_I_12 self.dry_I_13 = self.I_13_without_motor + self.motor.dry_I_13 @@ -696,18 +701,14 @@ def evaluate_inertias(self): CM_to_CPM = self.center_of_mass - self.center_of_propellant_position # Compute inertias - self.I_11 = ( - self.dry_I_11 - + self.motor.I_11 - + dry_mass * CM_to_CDM**2 - + prop_mass * CM_to_CPM**2 - ) - self.I_22 = ( - self.dry_I_22 - + self.motor.I_22 - + dry_mass * CM_to_CDM**2 - + prop_mass * CM_to_CPM**2 - ) + self.I_11 = parallel_axis_theorem_from_com( + self.dry_I_11, dry_mass, CM_to_CDM + ) + parallel_axis_theorem_from_com(self.motor.I_11, prop_mass, CM_to_CPM) + + self.I_22 = parallel_axis_theorem_from_com( + self.dry_I_22, dry_mass, CM_to_CDM + ) + parallel_axis_theorem_from_com(self.motor.I_22, prop_mass, CM_to_CPM) + self.I_33 = self.dry_I_33 + self.motor.I_33 self.I_12 = self.dry_I_12 + self.motor.I_12 self.I_13 = self.dry_I_13 + self.motor.I_13 diff --git a/rocketpy/tools.py b/rocketpy/tools.py index bcff91fb8..cb3094343 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -153,7 +153,7 @@ def time_num_to_date_string(time_num, units, timezone, calendar="gregorian"): """Convert time number (usually hours before a certain date) into two strings: one for the date (example: 2022.04.31) and one for the hour (example: 14). See cftime.num2date for details on units and calendar. - Automatically converts time number from UTC to local timezone based on + Automatically converts time number from UTC to local time zone based on lat, lon coordinates. This function was created originally for the EnvironmentAnalysis class. @@ -382,6 +382,33 @@ def check_requirement_version(module_name, version): return True +def parallel_axis_theorem_from_com(com_inertia_moment, mass, distance): + """Calculates the moment of inertia of a object relative to a new axis using + the parallel axis theorem. The new axis is parallel to and at a distance + 'distance' from the original axis, which *must* passes through the object's + center of mass. + + Parameters + ---------- + com_inertia_moment : float + Moment of inertia relative to the center of mass of the object. + mass : float + Mass of the object. + distance : float + Perpendicular distance between the original and new axis. + + Returns + ------- + float + Moment of inertia relative to the new axis. + + Reference + --------- + https://en.wikipedia.org/wiki/Parallel_axis_theorem + """ + return com_inertia_moment + mass * distance**2 + + # Flight