From c57c20072a6000605e549964becee03457541f26 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:09:26 -0300 Subject: [PATCH 01/25] ENH: add multiple argument support for shepard opt. - Minor bug in computation if distance zero solved --- rocketpy/mathutils/function.py | 80 +++++++++++++++++++--------------- 1 file changed, 46 insertions(+), 34 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 251d58199..8d97565ef 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -482,25 +482,33 @@ def get_value_opt(x): return y elif self.__interpolation__ == "shepard": - x_data = self.source[:, 0:-1] # Support for N-Dimensions - y_data = self.source[:, -1] - len_y_data = len(y_data) # A little speed up - # change the function's name to avoid mypy's error def get_value_opt_multiple(*args): - x = np.array([[float(val) for val in args]]) + x_data = self.source[:, 0:-1] # Support for N-Dimensions + y_data = self.source[:, -1] + + arg_stack = np.column_stack(args) + arg_qty, arg_dim = arg_stack.shape + result = np.zeros(arg_qty) + + # Reshape to vectorize calculations + x = arg_stack.reshape(arg_qty, 1, arg_dim) + sub_matrix = x_data - x - distances_squared = np.sum(sub_matrix**2, axis=1) + distances_squared = np.sum(sub_matrix**2, axis=2) - zero_distance_index = np.where(distances_squared == 0)[0] - if len(zero_distance_index) > 0: - return y_data[zero_distance_index[0]] + # Remove zero distances from further calculations + zero_distances = np.where(distances_squared == 0) + valid_indexes = np.ones(arg_qty, dtype=bool) + valid_indexes[zero_distances[0]] = False - weights = distances_squared ** (-1.5) - numerator_sum = np.sum(y_data * weights) - denominator_sum = np.sum(weights) + weights = distances_squared[valid_indexes] ** (-1.5) + numerator_sum = np.sum(y_data * weights, axis=1) + denominator_sum = np.sum(weights, axis=1) + result[valid_indexes] = numerator_sum / denominator_sum + result[~valid_indexes] = y_data[zero_distances[1]] - return numerator_sum / denominator_sum + return result get_value_opt = get_value_opt_multiple @@ -872,28 +880,32 @@ def get_value(self, *args): # Returns value for shepard interpolation elif self.__interpolation__ == "shepard": - if all(isinstance(arg, Iterable) for arg in args): - x = list(np.column_stack(args)) - else: - x = [[float(x) for x in list(args)]] - ans = x - x_data = self.source[:, 0:-1] + x_data = self.source[:, 0:-1] # Support for N-Dimensions y_data = self.source[:, -1] - for i, _ in enumerate(x): - numerator_sum = 0 - denominator_sum = 0 - for o, _ in enumerate(y_data): - sub = x_data[o] - x[i] - distance = (sub.dot(sub)) ** (0.5) - if distance == 0: - numerator_sum = y_data[o] - denominator_sum = 1 - break - weight = distance ** (-3) - numerator_sum = numerator_sum + y_data[o] * weight - denominator_sum = denominator_sum + weight - ans[i] = numerator_sum / denominator_sum - return ans if len(ans) > 1 else ans[0] + + arg_stack = np.column_stack(args) + arg_qty, arg_dim = arg_stack.shape + result = np.zeros(arg_qty) + + # Reshape to vectorize calculations + x = arg_stack.reshape(arg_qty, 1, arg_dim) + + sub_matrix = x_data - x + distances_squared = np.sum(sub_matrix**2, axis=2) + + # Remove zero distances from further calculations + zero_distances = np.where(distances_squared == 0) + valid_indexes = np.ones(arg_qty, dtype=bool) + valid_indexes[zero_distances[0]] = False + + weights = distances_squared[valid_indexes] ** (-1.5) + numerator_sum = np.sum(y_data * weights, axis=1) + denominator_sum = np.sum(weights, axis=1) + result[valid_indexes] = numerator_sum / denominator_sum + result[~valid_indexes] = y_data[zero_distances[1]] + + return result if len(result) > 1 else result[0] + # Returns value for polynomial interpolation function type elif self.__interpolation__ == "polynomial": if isinstance(args[0], (int, float)): From 3cd456f945c993afde829ff6c7d3802a6c2d8d08 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:10:01 -0300 Subject: [PATCH 02/25] TST: improve testing of higher dimensional functions. --- tests/test_function.py | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/tests/test_function.py b/tests/test_function.py index 183888b96..5b7daab5a 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -262,6 +262,7 @@ def test_multivariable_function_plot(mock_show): @pytest.mark.parametrize( "x,y,z_expected", [ + (0, 0, 1), (1, 0, 0), (0, 1, 0), (0, 0, 1), @@ -270,15 +271,40 @@ def test_multivariable_function_plot(mock_show): ([0, 0.5], [0, 0.5], [1, 1 / 3]), ], ) -def test_shepard_interpolation(x, y, z_expected): +def test_2d_shepard_interpolation(x, y, z_expected): """Test the shepard interpolation method of the Function class.""" # Test plane x + y + z = 1 source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] func = Function(source=source, inputs=["x", "y"], outputs=["z"]) z = func(x, y) + z_opt = func.get_value_opt(x, y) + assert np.isclose(z, z_opt, atol=1e-8).all() assert np.isclose(z, z_expected, atol=1e-8).all() +@pytest.mark.parametrize( + "x,y,z,w_expected", + [ + (0, 0, 0, 1), + (1, 0, 0, 0), + (0, 1, 0, 0), + (0, 0, 1, 0), + (0.5, 0.5, 0.5, 1 / 4), + (0.25, 0.25, 0.25, 0.700632626832), + ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), + ], +) +def test_3d_shepard_interpolation(x, y, z, w_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z + w = 1 + source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] + func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) + w = func(x, y, z) + w_opt = func.get_value_opt(x, y, z) + assert np.isclose(w, w_opt, atol=1e-8).all() + assert np.isclose(w_expected, w, atol=1e-8).all() + + @pytest.mark.parametrize("other", [1, 0.1, np.int_(1), np.float_(0.1), np.array([1])]) def test_sum_arithmetic_priority(other): """Test the arithmetic priority of the add operation of the Function class, From 954be1dad6c1dba92fa53ddf1fd8844cd750f549 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:22:51 -0300 Subject: [PATCH 03/25] FIX: interpolation return type if single value array. --- rocketpy/mathutils/function.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 8d97565ef..0e8213ac1 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -508,7 +508,7 @@ def get_value_opt_multiple(*args): result[valid_indexes] = numerator_sum / denominator_sum result[~valid_indexes] = y_data[zero_distances[1]] - return result + return result if len(result) > 1 else result[0] get_value_opt = get_value_opt_multiple From 7eab8ed751f1443070fbb410c99c38b75e7be90b Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:24:56 -0300 Subject: [PATCH 04/25] TST: move testing to unit module. --- tests/test_function.py | 46 ------------------------------------- tests/unit/test_function.py | 46 +++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 46 deletions(-) diff --git a/tests/test_function.py b/tests/test_function.py index 5b7daab5a..0b31d21ca 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -259,52 +259,6 @@ def test_multivariable_function_plot(mock_show): assert func.plot() == None -@pytest.mark.parametrize( - "x,y,z_expected", - [ - (0, 0, 1), - (1, 0, 0), - (0, 1, 0), - (0, 0, 1), - (0.5, 0.5, 1 / 3), - (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), - ([0, 0.5], [0, 0.5], [1, 1 / 3]), - ], -) -def test_2d_shepard_interpolation(x, y, z_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z = 1 - source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] - func = Function(source=source, inputs=["x", "y"], outputs=["z"]) - z = func(x, y) - z_opt = func.get_value_opt(x, y) - assert np.isclose(z, z_opt, atol=1e-8).all() - assert np.isclose(z, z_expected, atol=1e-8).all() - - -@pytest.mark.parametrize( - "x,y,z,w_expected", - [ - (0, 0, 0, 1), - (1, 0, 0, 0), - (0, 1, 0, 0), - (0, 0, 1, 0), - (0.5, 0.5, 0.5, 1 / 4), - (0.25, 0.25, 0.25, 0.700632626832), - ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), - ], -) -def test_3d_shepard_interpolation(x, y, z, w_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z + w = 1 - source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] - func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) - w = func(x, y, z) - w_opt = func.get_value_opt(x, y, z) - assert np.isclose(w, w_opt, atol=1e-8).all() - assert np.isclose(w_expected, w, atol=1e-8).all() - - @pytest.mark.parametrize("other", [1, 0.1, np.int_(1), np.float_(0.1), np.array([1])]) def test_sum_arithmetic_priority(other): """Test the arithmetic priority of the add operation of the Function class, diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index e41d30f04..ea52a88e6 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -48,6 +48,52 @@ def test_integral_spline_interpolation(request, func, a, b): ) +@pytest.mark.parametrize( + "x,y,z_expected", + [ + (0, 0, 1), + (1, 0, 0), + (0, 1, 0), + (0, 0, 1), + (0.5, 0.5, 1 / 3), + (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), + ([0, 0.5], [0, 0.5], [1, 1 / 3]), + ], +) +def test_2d_shepard_interpolation(x, y, z_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z = 1 + source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] + func = Function(source=source, inputs=["x", "y"], outputs=["z"]) + z = func(x, y) + z_opt = func.get_value_opt(x, y) + assert np.isclose(z, z_opt, atol=1e-8).all() + assert np.isclose(z, z_expected, atol=1e-8).all() + + +@pytest.mark.parametrize( + "x,y,z,w_expected", + [ + (0, 0, 0, 1), + (1, 0, 0, 0), + (0, 1, 0, 0), + (0, 0, 1, 0), + (0.5, 0.5, 0.5, 1 / 4), + (0.25, 0.25, 0.25, 0.700632626832), + ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), + ], +) +def test_3d_shepard_interpolation(x, y, z, w_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z + w = 1 + source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] + func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) + w = func(x, y, z) + w_opt = func.get_value_opt(x, y, z) + assert np.isclose(w, w_opt, atol=1e-8).all() + assert np.isclose(w_expected, w, atol=1e-8).all() + + @pytest.mark.parametrize( "func_input, derivative_input, expected_derivative", [ From 8fdfa555430795895be5eedd6710f70468928ba1 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Mon, 18 Dec 2023 21:29:06 -0300 Subject: [PATCH 05/25] MNT: update CHANGELOG with changes. --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a8cbb9ed6..33926e5ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,7 +32,7 @@ straightforward as possible. ### Added -- +- ENH: Shepard Optimized Interpolation - Multiple Inputs Support [#515](https://github.com/RocketPy-Team/RocketPy/pull/515) ### Changed From 34ed384b50aff46a02173d98d8f3d34024d20769 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Tue, 19 Dec 2023 14:59:53 -0300 Subject: [PATCH 06/25] MNT: solve minor test placement issue. --- tests/test_function.py | 50 +++++++++++++++++++++++++++++++++++++ tests/unit/test_function.py | 46 ---------------------------------- 2 files changed, 50 insertions(+), 46 deletions(-) diff --git a/tests/test_function.py b/tests/test_function.py index 0b31d21ca..9cc2e68f0 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -215,6 +215,56 @@ def test_multivariable_dataset(a, b): assert np.isclose(func(a, b), a + b, atol=1e-6) +@pytest.mark.parametrize( + "x,y,z_expected", + [ + (0, 0, 1), + (1, 0, 0), + (0, 1, 0), + (0, 0, 1), + (0.5, 0.5, 1 / 3), + (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), + ([0, 0.5], [0, 0.5], [1, 1 / 3]), + ], +) +def test_2d_shepard_interpolation(x, y, z_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z = 1 + source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] + func = Function( + source=source, inputs=["x", "y"], outputs=["z"], interpolation="shepard" + ) + z = func(x, y) + z_opt = func.get_value_opt(x, y) + assert np.isclose(z, z_opt, atol=1e-8).all() + assert np.isclose(z_expected, z, atol=1e-8).all() + + +@pytest.mark.parametrize( + "x,y,z,w_expected", + [ + (0, 0, 0, 1), + (1, 0, 0, 0), + (0, 1, 0, 0), + (0, 0, 1, 0), + (0.5, 0.5, 0.5, 1 / 4), + (0.25, 0.25, 0.25, 0.700632626832), + ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), + ], +) +def test_3d_shepard_interpolation(x, y, z, w_expected): + """Test the shepard interpolation method of the Function class.""" + # Test plane x + y + z + w = 1 + source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] + func = Function( + source=source, inputs=["x", "y", "z"], outputs=["w"], interpolation="shepard" + ) + w = func(x, y, z) + w_opt = func.get_value_opt(x, y, z) + assert np.isclose(w, w_opt, atol=1e-8).all() + assert np.isclose(w_expected, w, atol=1e-8).all() + + @pytest.mark.parametrize("a", [-1, -0.5, 0, 0.5, 1]) @pytest.mark.parametrize("b", [-1, -0.5, 0, 0.5, 1]) def test_multivariable_function(a, b): diff --git a/tests/unit/test_function.py b/tests/unit/test_function.py index ea52a88e6..e41d30f04 100644 --- a/tests/unit/test_function.py +++ b/tests/unit/test_function.py @@ -48,52 +48,6 @@ def test_integral_spline_interpolation(request, func, a, b): ) -@pytest.mark.parametrize( - "x,y,z_expected", - [ - (0, 0, 1), - (1, 0, 0), - (0, 1, 0), - (0, 0, 1), - (0.5, 0.5, 1 / 3), - (0.25, 0.25, 25 / (25 + 2 * 5**0.5)), - ([0, 0.5], [0, 0.5], [1, 1 / 3]), - ], -) -def test_2d_shepard_interpolation(x, y, z_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z = 1 - source = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] - func = Function(source=source, inputs=["x", "y"], outputs=["z"]) - z = func(x, y) - z_opt = func.get_value_opt(x, y) - assert np.isclose(z, z_opt, atol=1e-8).all() - assert np.isclose(z, z_expected, atol=1e-8).all() - - -@pytest.mark.parametrize( - "x,y,z,w_expected", - [ - (0, 0, 0, 1), - (1, 0, 0, 0), - (0, 1, 0, 0), - (0, 0, 1, 0), - (0.5, 0.5, 0.5, 1 / 4), - (0.25, 0.25, 0.25, 0.700632626832), - ([0, 0.5], [0, 0.5], [0, 0.5], [1, 1 / 4]), - ], -) -def test_3d_shepard_interpolation(x, y, z, w_expected): - """Test the shepard interpolation method of the Function class.""" - # Test plane x + y + z + w = 1 - source = [(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)] - func = Function(source=source, inputs=["x", "y", "z"], outputs=["w"]) - w = func(x, y, z) - w_opt = func.get_value_opt(x, y, z) - assert np.isclose(w, w_opt, atol=1e-8).all() - assert np.isclose(w_expected, w, atol=1e-8).all() - - @pytest.mark.parametrize( "func_input, derivative_input, expected_derivative", [ From 4ed16bee1d88fe781b6fad5bff73aa9d49720898 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Tue, 19 Dec 2023 18:41:00 -0300 Subject: [PATCH 07/25] MNT: avoid code interpolation code repetition. --- rocketpy/mathutils/function.py | 93 ++++++++++++++++------------------ 1 file changed, 43 insertions(+), 50 deletions(-) diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 0e8213ac1..ceeb6a2e2 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -484,31 +484,7 @@ def get_value_opt(x): elif self.__interpolation__ == "shepard": # change the function's name to avoid mypy's error def get_value_opt_multiple(*args): - x_data = self.source[:, 0:-1] # Support for N-Dimensions - y_data = self.source[:, -1] - - arg_stack = np.column_stack(args) - arg_qty, arg_dim = arg_stack.shape - result = np.zeros(arg_qty) - - # Reshape to vectorize calculations - x = arg_stack.reshape(arg_qty, 1, arg_dim) - - sub_matrix = x_data - x - distances_squared = np.sum(sub_matrix**2, axis=2) - - # Remove zero distances from further calculations - zero_distances = np.where(distances_squared == 0) - valid_indexes = np.ones(arg_qty, dtype=bool) - valid_indexes[zero_distances[0]] = False - - weights = distances_squared[valid_indexes] ** (-1.5) - numerator_sum = np.sum(y_data * weights, axis=1) - denominator_sum = np.sum(weights, axis=1) - result[valid_indexes] = numerator_sum / denominator_sum - result[~valid_indexes] = y_data[zero_distances[1]] - - return result if len(result) > 1 else result[0] + return self.__interpolate_shepard__(args) get_value_opt = get_value_opt_multiple @@ -880,31 +856,7 @@ def get_value(self, *args): # Returns value for shepard interpolation elif self.__interpolation__ == "shepard": - x_data = self.source[:, 0:-1] # Support for N-Dimensions - y_data = self.source[:, -1] - - arg_stack = np.column_stack(args) - arg_qty, arg_dim = arg_stack.shape - result = np.zeros(arg_qty) - - # Reshape to vectorize calculations - x = arg_stack.reshape(arg_qty, 1, arg_dim) - - sub_matrix = x_data - x - distances_squared = np.sum(sub_matrix**2, axis=2) - - # Remove zero distances from further calculations - zero_distances = np.where(distances_squared == 0) - valid_indexes = np.ones(arg_qty, dtype=bool) - valid_indexes[zero_distances[0]] = False - - weights = distances_squared[valid_indexes] ** (-1.5) - numerator_sum = np.sum(y_data * weights, axis=1) - denominator_sum = np.sum(weights, axis=1) - result[valid_indexes] = numerator_sum / denominator_sum - result[~valid_indexes] = y_data[zero_distances[1]] - - return result if len(result) > 1 else result[0] + return self.__interpolate_shepard__(args) # Returns value for polynomial interpolation function type elif self.__interpolation__ == "polynomial": @@ -1656,6 +1608,47 @@ def __interpolate_akima__(self): coeffs[4 * i : 4 * i + 4] = np.linalg.solve(matrix, result) self.__akima_coefficients__ = coeffs + def __interpolate_shepard__(self, args): + """Calculates the shepard interpolation from the given arguments. + The shepard interpolation is computed by a inverse distance weighting + in a vectorized manner. + + Parameters + ---------- + args : scalar, list + Values where the Function is to be evaluated. + + Returns + ------- + result : scalar, list + The result of the interpolation. + """ + x_data = self.source[:, 0:-1] # Support for N-Dimensions + y_data = self.source[:, -1] + + arg_stack = np.column_stack(args) + arg_qty, arg_dim = arg_stack.shape + result = np.zeros(arg_qty) + + # Reshape to vectorize calculations + x = arg_stack.reshape(arg_qty, 1, arg_dim) + + sub_matrix = x_data - x + distances_squared = np.sum(sub_matrix**2, axis=2) + + # Remove zero distances from further calculations + zero_distances = np.where(distances_squared == 0) + valid_indexes = np.ones(arg_qty, dtype=bool) + valid_indexes[zero_distances[0]] = False + + weights = distances_squared[valid_indexes] ** (-1.5) + numerator_sum = np.sum(y_data * weights, axis=1) + denominator_sum = np.sum(weights, axis=1) + result[valid_indexes] = numerator_sum / denominator_sum + result[~valid_indexes] = y_data[zero_distances[1]] + + return result if len(result) > 1 else result[0] + def __neg__(self): """Negates the Function object. The result has the same effect as multiplying the Function by -1. From acbc14dae924d874e6a2803e2f494be92eec7f09 Mon Sep 17 00:00:00 2001 From: Pedro Bressan Date: Thu, 18 Jan 2024 20:56:04 -0300 Subject: [PATCH 08/25] MNT: remove repeated test case. --- tests/test_function.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_function.py b/tests/test_function.py index 9cc2e68f0..352feb7b3 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -218,7 +218,6 @@ def test_multivariable_dataset(a, b): @pytest.mark.parametrize( "x,y,z_expected", [ - (0, 0, 1), (1, 0, 0), (0, 1, 0), (0, 0, 1), From 5c1267f4db57198b997f1b993c2df0de332dc954 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 15 Dec 2023 20:16:18 +0100 Subject: [PATCH 09/25] ENH: add barometric_height to environment --- rocketpy/environment/environment.py | 56 ++++++++++++++++++++++++++++- 1 file changed, 55 insertions(+), 1 deletion(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 63d24fa9e..adf030806 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -111,6 +111,10 @@ class Environment: Environment.pressure : Function Air pressure in Pa as a function of altitude. Can be accessed as regular array, or called as a Function. See Function for more information. + Environment.barometric_height : Function + Geometric height above sea level in m as a function of pressure. Can be + accessed as regular array, or called as a Function. See Function for + more information. Environment.temperature : Function Air temperature in K as a function of altitude. Can be accessed as regular array, or called as a Function. See Function for more @@ -1315,6 +1319,8 @@ def process_standard_atmosphere(self): # Save temperature, pressure and wind profiles self.pressure = self.pressure_ISA + self.barometric_height = self.barometric_height_ISA + self.temperature = self.temperature_ISA self.wind_direction = Function( 0, @@ -1433,6 +1439,7 @@ def process_custom_atmosphere( if pressure is None: # Use standard atmosphere self.pressure = self.pressure_ISA + self.barometric_height = self.barometric_height_ISA else: # Use custom input self.pressure = Function( @@ -1441,6 +1448,11 @@ def process_custom_atmosphere( outputs="Pressure (Pa)", interpolation="linear", ) + self.barometric_height = self.pressure.inverse_function().set_discrete( + 0, max_expected_height, 1000, extrapolation="constant" + ) + self.barometric_height.set_inputs("Pressure (Pa)") + self.barometric_height.set_outputs("Height Above Sea Level (m)") # Check maximum height of custom pressure input if not callable(self.pressure.source): max_expected_height = max(self.pressure[-1, 0], max_expected_height) @@ -1605,6 +1617,12 @@ def process_windy_atmosphere(self, model="ECMWF"): outputs="Pressure (Pa)", interpolation="linear", ) + self.barometric_height = Function( + data_array[:, (0, 1)], + inputs="Pressure (Pa)", + outputs="Height Above Sea Level (m)", + interpolation="linear", + ) self.temperature = Function( data_array[:, (1, 2)], inputs="Height Above Sea Level (m)", @@ -1732,6 +1750,12 @@ def process_wyoming_sounding(self, file): outputs="Pressure (Pa)", interpolation="linear", ) + self.barometric_height = Function( + data_array[:, (0, 1)], + inputs="Pressure (Pa)", + outputs="Height Above Sea Level (m)", + interpolation="linear", + ) # Retrieve temperature from data array data_array[:, 2] = data_array[:, 2] + 273.15 # Converts C to K @@ -1845,6 +1869,7 @@ def process_noaaruc_sounding(self, file): # Extract pressure as a function of height pressure_array = [] + barometric_height_array = [] for line in lines: # Split line into columns columns = re.split(" +", line)[1:] @@ -1858,7 +1883,9 @@ def process_noaaruc_sounding(self, file): if max(columns) != 99999: # Save value pressure_array.append(columns) + barometric_height_array.append([columns[1], columns[0]]) pressure_array = np.array(pressure_array) + barometric_height_array = np.array(barometric_height_array) # Extract temperature as a function of height temperature_array = [] @@ -1905,6 +1932,14 @@ def process_noaaruc_sounding(self, file): outputs="Pressure (Pa)", interpolation="linear", ) + # construct barometric height array from pressure array + barometric_height_array[:, 0] = 10 * barometric_height_array[:, 0] + self.barometric_height = Function( + barometric_height_array, + inputs="Pressure (Pa)", + outputs="Height Above Sea Level (m)", + interpolation="linear", + ) # Convert 10*C to K and save values temperature_array[:, 1] = ( @@ -2274,6 +2309,12 @@ def process_forecast_reanalysis(self, file, dictionary): outputs="Pressure (Pa)", interpolation="linear", ) + self.barometric_height = Function( + data_array[:, (0, 1)], + inputs="Pressure (Pa)", + outputs="Height Above Sea Level (m)", + interpolation="linear", + ) self.temperature = Function( data_array[:, (1, 2)], inputs="Height Above Sea Level (m)", @@ -2803,6 +2844,12 @@ def select_ensemble_member(self, member=0): outputs="Pressure (Pa)", interpolation="linear", ) + self.barometric_height = Function( + data_array[:, (0, 1)], + inputs="Pressure (Pa)", + outputs="Height Above Sea Level (m)", + interpolation="linear", + ) self.temperature = Function( data_array[:, (1, 2)], inputs="Height Above Sea Level (m)", @@ -2965,7 +3012,14 @@ def pressure_function(h): outputs="Pressure (Pa)", ) - return None + # Save international standard atmosphere height by pressure profile + # and discretize it. This is done to speed up the calculations in the + # trajectory simulation. + self.barometric_height_ISA = self.pressure_ISA.inverse_function().set_discrete( + pressure[-1], pressure[0], 1000, extrapolation="constant" + ) + self.barometric_height_ISA.set_inputs("Pressure (Pa)") + self.barometric_height_ISA.set_outputs("Height Above Sea Level (m)") def calculate_density_profile(self): """Compute the density of the atmosphere as a function of From f301352528c25edeb5137e474fe4d70b3f8186d7 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 15 Dec 2023 20:24:26 +0100 Subject: [PATCH 10/25] ENH: get hAGL from barometric_height --- rocketpy/simulation/flight.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index c4a719e88..299ddaead 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -697,9 +697,8 @@ def __init__( parachute.noisy_pressure_signal.append([node.t, pressure + noise]) # Gets height above ground level considering noise hAGL = ( - self.env.pressure.find_input( + self.env.barometric_height( pressure + noise, - self.y_sol[2], ) - self.env.elevation ) @@ -1006,9 +1005,8 @@ def __init__( ) # Gets height above ground level considering noise hAGL = ( - self.env.pressure.find_input( + self.env.barometric_height( pressure + noise, - overshootable_node.y[2], ) - self.env.elevation ) From ead932275e386997276cc8d1ccd7d97d09dd4e61 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Fri, 15 Dec 2023 20:39:55 +0100 Subject: [PATCH 11/25] MNT: update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index d4cd4be85..4c5ad55ad 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -39,6 +39,7 @@ straightforward as possible. - MNT: Add repr method to Parachute class [#490](https://github.com/RocketPy-Team/RocketPy/pull/490) - ENH: Function Reverse Arithmetic Priority [#488](https://github.com/RocketPy-Team/RocketPy/pull/488) - DOC: Update Header Related Docs +- ENH Precalculate Barometric Height [#511](https://github.com/RocketPy-Team/RocketPy/pull/511) ### Fixed From 07d4c338ba6b5b0fcdbdac3bda81f87b0d92b136 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 18 Dec 2023 04:11:29 +0100 Subject: [PATCH 12/25] DOC: improve black formatting --- rocketpy/simulation/flight.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 299ddaead..85f0edc25 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -697,9 +697,7 @@ def __init__( parachute.noisy_pressure_signal.append([node.t, pressure + noise]) # Gets height above ground level considering noise hAGL = ( - self.env.barometric_height( - pressure + noise, - ) + self.env.barometric_height(pressure + noise) - self.env.elevation ) if parachute.triggerfunc(pressure + noise, hAGL, self.y_sol): @@ -1005,9 +1003,7 @@ def __init__( ) # Gets height above ground level considering noise hAGL = ( - self.env.barometric_height( - pressure + noise, - ) + self.env.barometric_height(pressure + noise) - self.env.elevation ) From 2b06a7e42f0353deb3a53c2e9ffb1b1d8f66c926 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 18 Dec 2023 19:26:50 +0100 Subject: [PATCH 13/25] ENH: add extrapolation to barometric_height --- rocketpy/environment/environment.py | 88 ++++++++++++++++++++++++++--- 1 file changed, 79 insertions(+), 9 deletions(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index adf030806..94857a961 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -1449,7 +1449,7 @@ def process_custom_atmosphere( interpolation="linear", ) self.barometric_height = self.pressure.inverse_function().set_discrete( - 0, max_expected_height, 1000, extrapolation="constant" + 0, max_expected_height, 100, extrapolation="constant" ) self.barometric_height.set_inputs("Pressure (Pa)") self.barometric_height.set_outputs("Height Above Sea Level (m)") @@ -1617,8 +1617,21 @@ def process_windy_atmosphere(self, model="ECMWF"): outputs="Pressure (Pa)", interpolation="linear", ) + # Linearly extrapolate pressure to ground level + columns = data_array[:, (0, 1)] + pressure_ground = np.array( + [ + [ + columns[1][0] + + (-columns[1][1] / (columns[0][1] - columns[1][1])) + * (columns[0][0] - columns[1][0]), + 0, + ] + ] + ) + bar_height = np.concatenate((pressure_ground, columns)) self.barometric_height = Function( - data_array[:, (0, 1)], + bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", @@ -1750,8 +1763,21 @@ def process_wyoming_sounding(self, file): outputs="Pressure (Pa)", interpolation="linear", ) + # Linearly extrapolate pressure to ground level + columns = data_array[:, (0, 1)] + pressure_ground = np.array( + [ + [ + columns[1][0] + + (-columns[1][1] / (columns[0][1] - columns[1][1])) + * (columns[0][0] - columns[1][0]), + 0, + ] + ] + ) + bar_height = np.concatenate((pressure_ground, columns)) self.barometric_height = Function( - data_array[:, (0, 1)], + bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", @@ -1932,10 +1958,28 @@ def process_noaaruc_sounding(self, file): outputs="Pressure (Pa)", interpolation="linear", ) - # construct barometric height array from pressure array + # Converts 10*hPa to Pa and save values barometric_height_array[:, 0] = 10 * barometric_height_array[:, 0] + # Linearly extrapolate pressure to ground level + pressure_ground = np.array( + [ + [ + barometric_height_array[1][0] + + ( + -barometric_height_array[1][1] + / ( + barometric_height_array[0][1] + - barometric_height_array[1][1] + ) + ) + * (barometric_height_array[0][0] - barometric_height_array[1][0]), + 0, + ] + ] + ) + bar_height = np.concatenate((pressure_ground, barometric_height_array)) self.barometric_height = Function( - barometric_height_array, + bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", @@ -2309,8 +2353,21 @@ def process_forecast_reanalysis(self, file, dictionary): outputs="Pressure (Pa)", interpolation="linear", ) + # Linearly extrapolate pressure to ground level + columns = data_array[:, (0, 1)] + pressure_ground = np.array( + [ + [ + columns[1][0] + + (-columns[1][1] / (columns[0][1] - columns[1][1])) + * (columns[0][0] - columns[1][0]), + 0, + ] + ] + ) + bar_height = np.concatenate((pressure_ground, columns)) self.barometric_height = Function( - data_array[:, (0, 1)], + bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", @@ -2844,11 +2901,24 @@ def select_ensemble_member(self, member=0): outputs="Pressure (Pa)", interpolation="linear", ) + # Linearly extrapolate pressure to ground level + columns = data_array[:, (0, 1)] + pressure_ground = np.array( + [ + [ + columns[1][0] + + (-columns[1][1] / (columns[0][1] - columns[1][1])) + * (columns[0][0] - columns[1][0]), + 0, + ] + ] + ) + bar_height = np.concatenate((pressure_ground, columns)) self.barometric_height = Function( - data_array[:, (0, 1)], + bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", - interpolation="linear", + interpolation="linearS", ) self.temperature = Function( data_array[:, (1, 2)], @@ -3016,7 +3086,7 @@ def pressure_function(h): # and discretize it. This is done to speed up the calculations in the # trajectory simulation. self.barometric_height_ISA = self.pressure_ISA.inverse_function().set_discrete( - pressure[-1], pressure[0], 1000, extrapolation="constant" + pressure[-1], pressure[0], 100, extrapolation="constant" ) self.barometric_height_ISA.set_inputs("Pressure (Pa)") self.barometric_height_ISA.set_outputs("Height Above Sea Level (m)") From e787b1180c15bba8cd21d8b78f074937d271a76d Mon Sep 17 00:00:00 2001 From: MateusStano Date: Mon, 18 Dec 2023 19:27:03 +0100 Subject: [PATCH 14/25] TST: add barometric_height to tests --- tests/test_environment.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_environment.py b/tests/test_environment.py index 7a23828c2..db80eea6d 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -20,6 +20,7 @@ def test_standard_atmosphere(mock_show, example_env): assert example_env.info() == None assert example_env.all_info() == None assert abs(example_env.pressure(0) - 101325.0) < 1e-8 + assert example_env.barometric_height(101325.0) < 1e-2 assert example_env.prints.print_earth_details() == None @@ -43,6 +44,7 @@ def test_custom_atmosphere(mock_show, example_env): ) assert example_env.all_info() == None assert abs(example_env.pressure(0) - 101325.0) < 1e-8 + assert example_env.barometric_height(101325.0) < 1e-2 assert abs(example_env.wind_velocity_x(0) - 5) < 1e-8 assert abs(example_env.temperature(100) - 300) < 1e-8 @@ -71,6 +73,7 @@ def test_wyoming_sounding_atmosphere(mock_show, example_env): pass assert example_env.all_info() == None assert abs(example_env.pressure(0) - 93600.0) < 1e-8 + assert example_env.barometric_height(example_env.pressure(0)) == 722.0 assert abs(example_env.wind_velocity_x(0) - -2.9005178894925043) < 1e-8 assert abs(example_env.temperature(100) - 291.75) < 1e-8 From 0551317841b232f224876e84ce383f3f02a71962 Mon Sep 17 00:00:00 2001 From: MateusStano Date: Thu, 18 Jan 2024 19:12:11 +0100 Subject: [PATCH 15/25] ENH: add use of natural extrapolation to barometric_height --- rocketpy/environment/environment.py | 79 ++++------------------------- 1 file changed, 11 insertions(+), 68 deletions(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 94857a961..2778f19fb 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -1618,23 +1618,13 @@ def process_windy_atmosphere(self, model="ECMWF"): interpolation="linear", ) # Linearly extrapolate pressure to ground level - columns = data_array[:, (0, 1)] - pressure_ground = np.array( - [ - [ - columns[1][0] - + (-columns[1][1] / (columns[0][1] - columns[1][1])) - * (columns[0][0] - columns[1][0]), - 0, - ] - ] - ) - bar_height = np.concatenate((pressure_ground, columns)) + bar_height = data_array[:, (0, 1)] self.barometric_height = Function( bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", + extrapolation="natural", ) self.temperature = Function( data_array[:, (1, 2)], @@ -1764,23 +1754,13 @@ def process_wyoming_sounding(self, file): interpolation="linear", ) # Linearly extrapolate pressure to ground level - columns = data_array[:, (0, 1)] - pressure_ground = np.array( - [ - [ - columns[1][0] - + (-columns[1][1] / (columns[0][1] - columns[1][1])) - * (columns[0][0] - columns[1][0]), - 0, - ] - ] - ) - bar_height = np.concatenate((pressure_ground, columns)) + bar_height = data_array[:, (0, 1)] self.barometric_height = Function( bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", + extrapolation="natural", ) # Retrieve temperature from data array @@ -1960,29 +1940,12 @@ def process_noaaruc_sounding(self, file): ) # Converts 10*hPa to Pa and save values barometric_height_array[:, 0] = 10 * barometric_height_array[:, 0] - # Linearly extrapolate pressure to ground level - pressure_ground = np.array( - [ - [ - barometric_height_array[1][0] - + ( - -barometric_height_array[1][1] - / ( - barometric_height_array[0][1] - - barometric_height_array[1][1] - ) - ) - * (barometric_height_array[0][0] - barometric_height_array[1][0]), - 0, - ] - ] - ) - bar_height = np.concatenate((pressure_ground, barometric_height_array)) self.barometric_height = Function( - bar_height, + barometric_height_array, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", + extrapolation="natural", ) # Convert 10*C to K and save values @@ -2354,23 +2317,13 @@ def process_forecast_reanalysis(self, file, dictionary): interpolation="linear", ) # Linearly extrapolate pressure to ground level - columns = data_array[:, (0, 1)] - pressure_ground = np.array( - [ - [ - columns[1][0] - + (-columns[1][1] / (columns[0][1] - columns[1][1])) - * (columns[0][0] - columns[1][0]), - 0, - ] - ] - ) - bar_height = np.concatenate((pressure_ground, columns)) + bar_height = data_array[:, (0, 1)] self.barometric_height = Function( bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", interpolation="linear", + extrapolation="natural", ) self.temperature = Function( data_array[:, (1, 2)], @@ -2902,23 +2855,13 @@ def select_ensemble_member(self, member=0): interpolation="linear", ) # Linearly extrapolate pressure to ground level - columns = data_array[:, (0, 1)] - pressure_ground = np.array( - [ - [ - columns[1][0] - + (-columns[1][1] / (columns[0][1] - columns[1][1])) - * (columns[0][0] - columns[1][0]), - 0, - ] - ] - ) - bar_height = np.concatenate((pressure_ground, columns)) + bar_height = data_array[:, (0, 1)] self.barometric_height = Function( bar_height, inputs="Pressure (Pa)", outputs="Height Above Sea Level (m)", - interpolation="linearS", + interpolation="linear", + extrapolation="natural", ) self.temperature = Function( data_array[:, (1, 2)], From 4225405d6ae10f78608c95d6c0ec6ef467844d7e Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sat, 20 Jan 2024 00:10:40 -0500 Subject: [PATCH 16/25] MNT: Add new variable barometric_height_ISA to .pylintrc --- .pylintrc | 1 + 1 file changed, 1 insertion(+) diff --git a/.pylintrc b/.pylintrc index 97a3abeb1..d6116873f 100644 --- a/.pylintrc +++ b/.pylintrc @@ -187,6 +187,7 @@ function-naming-style=snake_case # Good variable names which should always be accepted, separated by a comma. good-names=FlightPhases, WindroseAxes, + barometric_height_ISA, # Good variable names regexes, separated by a comma. If names match any regex, From c1216b25f71a446c5c121fe79a027ec35613c61a Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Mon, 22 Jan 2024 15:50:41 -0500 Subject: [PATCH 17/25] MNT: Add quaternion conversion functions to rocketpy.tools and update Flight --- .pylintrc | 1 + CHANGELOG.md | 1 + rocketpy/simulation/flight.py | 36 +++++++++----------- rocketpy/tools.py | 64 +++++++++++++++++++++++++++++++++++ 4 files changed, 82 insertions(+), 20 deletions(-) diff --git a/.pylintrc b/.pylintrc index 2611847f2..8cebb954c 100644 --- a/.pylintrc +++ b/.pylintrc @@ -436,6 +436,7 @@ disable=raw-checker-failed, no-else-return, inconsistent-return-statements, unspecified-encoding, + no-member, # because we use funcify_method decorator # Enable the message, report, category or checker with the given id(s). You can # either give multiple identifier separated by comma (,) or put this option diff --git a/CHANGELOG.md b/CHANGELOG.md index ad4ac536b..a3f73a2fe 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,7 @@ straightforward as possible. - MNT: Add repr method to Parachute class [#490](https://github.com/RocketPy-Team/RocketPy/pull/490) - ENH: Function Reverse Arithmetic Priority [#488](https://github.com/RocketPy-Team/RocketPy/pull/488) - DOC: Update Header Related Docs +- MNT: Encapsulate quaternion conversions [#537](https://github.com/RocketPy-Team/RocketPy/pull/537) ### Fixed diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 1c2549dfa..0b6bf3c10 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -11,7 +11,12 @@ from ..mathutils.vector_matrix import Matrix, Vector from ..plots.flight_plots import _FlightPlots from ..prints.flight_prints import _FlightPrints -from ..tools import find_closest +from ..tools import ( + find_closest, + quaternions_to_phi, + quaternions_to_precession, + quaternions_to_theta, +) class Flight: @@ -2295,33 +2300,24 @@ def lateral_attitude_angle(self): @funcify_method("Time (s)", "Precession Angle - ψ (°)", "spline", "constant") def psi(self): """Precession angle as a Function of time.""" - psi = (180 / np.pi) * ( - np.arctan2(self.e3[:, 1], self.e0[:, 1]) - + np.arctan2(-self.e2[:, 1], -self.e1[:, 1]) - ) # Precession angle - psi = np.column_stack([self.time, psi]) # Precession angle - return psi + psi = quaternions_to_precession( + self.e0.y_array, self.e1.y_array, self.e2.y_array, self.e3.y_array + ) + return np.column_stack([self.time, psi]) @funcify_method("Time (s)", "Spin Angle - φ (°)", "spline", "constant") def phi(self): """Spin angle as a Function of time.""" - phi = (180 / np.pi) * ( - np.arctan2(self.e3[:, 1], self.e0[:, 1]) - - np.arctan2(-self.e2[:, 1], -self.e1[:, 1]) - ) # Spin angle - phi = np.column_stack([self.time, phi]) # Spin angle - return phi + phi = quaternions_to_phi( + self.e0.y_array, self.e1.y_array, self.e2.y_array, self.e3.y_array + ) + return np.column_stack([self.time, phi]) @funcify_method("Time (s)", "Nutation Angle - θ (°)", "spline", "constant") def theta(self): """Nutation angle as a Function of time.""" - theta = ( - (180 / np.pi) - * 2 - * np.arcsin(-((self.e1[:, 1] ** 2 + self.e2[:, 1] ** 2) ** 0.5)) - ) # Nutation angle - theta = np.column_stack([self.time, theta]) # Nutation angle - return theta + theta = quaternions_to_theta(self.e1.y_array, self.e2.y_array) + return np.column_stack([self.time, theta]) # Fluid Mechanics variables # Freestream Velocity diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 6e765ca7a..21ac0006a 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -3,6 +3,7 @@ import re from bisect import bisect_left +import numpy as np import pytz from cftime import num2pydate from packaging import version as packaging_version @@ -381,6 +382,69 @@ def check_requirement_version(module_name, version): return True +# Flight + + +def quaternions_to_precession(e0, e1, e2, e3): + """Calculates the Precession angle + + Parameters + ---------- + e0 : float + Euler parameter 0, must be between -1 and 1 + e1 : float + Euler parameter 1, must be between -1 and 1 + e2 : float + Euler parameter 2, must be between -1 and 1 + e3 : float + Euler parameter 3, must be between -1 and 1 + Returns + ------- + float + Euler Precession angle in degrees + """ + return (180 / np.pi) * (np.arctan2(e3, e0) + np.arctan2(-e2, -e1)) + + +def quaternions_to_phi(e0, e1, e2, e3): + """Calculates the Spin angle from quaternions. + + Parameters + ---------- + e0 : float + Euler parameter 0, must be between -1 and 1 + e1 : float + Euler parameter 1, must be between -1 and 1 + e2 : float + Euler parameter 2, must be between -1 and 1 + e3 : float + Euler parameter 3, must be between -1 and 1 + + Returns + ------- + float + Euler Spin angle in degrees + """ + return (180 / np.pi) * (np.arctan2(e3, e0) - np.arctan2(-e2, -e1)) + + +def quaternions_to_theta(e1, e2): + """Calculates the Nutation angle from quaternions. + + Parameters + ---------- + e1 : float + Euler parameter 1, must be between -1 and 1 + e2 : float + Euler parameter 2, must be between -1 and 1 + Returns + ------- + float + Euler Nutation angle in degrees + """ + return (180 / np.pi) * 2 * np.arcsin(-((e1**2 + e2**2) ** 0.5)) + + if __name__ == "__main__": import doctest From cdb42ece907a26f6b941eb404118607396214f76 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> Date: Tue, 23 Jan 2024 10:43:16 -0500 Subject: [PATCH 18/25] MNT: rename quaternions conversion functions --- rocketpy/simulation/flight.py | 8 ++++---- rocketpy/tools.py | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 0b6bf3c10..1d9017152 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -13,9 +13,9 @@ from ..prints.flight_prints import _FlightPrints from ..tools import ( find_closest, - quaternions_to_phi, + quaternions_to_spin, quaternions_to_precession, - quaternions_to_theta, + quaternions_to_nutation, ) @@ -2308,7 +2308,7 @@ def psi(self): @funcify_method("Time (s)", "Spin Angle - φ (°)", "spline", "constant") def phi(self): """Spin angle as a Function of time.""" - phi = quaternions_to_phi( + phi = quaternions_to_spin( self.e0.y_array, self.e1.y_array, self.e2.y_array, self.e3.y_array ) return np.column_stack([self.time, phi]) @@ -2316,7 +2316,7 @@ def phi(self): @funcify_method("Time (s)", "Nutation Angle - θ (°)", "spline", "constant") def theta(self): """Nutation angle as a Function of time.""" - theta = quaternions_to_theta(self.e1.y_array, self.e2.y_array) + theta = quaternions_to_nutation(self.e1.y_array, self.e2.y_array) return np.column_stack([self.time, theta]) # Fluid Mechanics variables diff --git a/rocketpy/tools.py b/rocketpy/tools.py index 21ac0006a..bcff91fb8 100644 --- a/rocketpy/tools.py +++ b/rocketpy/tools.py @@ -406,7 +406,7 @@ def quaternions_to_precession(e0, e1, e2, e3): return (180 / np.pi) * (np.arctan2(e3, e0) + np.arctan2(-e2, -e1)) -def quaternions_to_phi(e0, e1, e2, e3): +def quaternions_to_spin(e0, e1, e2, e3): """Calculates the Spin angle from quaternions. Parameters @@ -428,7 +428,7 @@ def quaternions_to_phi(e0, e1, e2, e3): return (180 / np.pi) * (np.arctan2(e3, e0) - np.arctan2(-e2, -e1)) -def quaternions_to_theta(e1, e2): +def quaternions_to_nutation(e1, e2): """Calculates the Nutation angle from quaternions. Parameters From d3e7e318cc76697d6d1af46d13355cdb44608b48 Mon Sep 17 00:00:00 2001 From: MateusStano <69485049+MateusStano@users.noreply.github.com> Date: Thu, 25 Jan 2024 12:35:38 -0300 Subject: [PATCH 19/25] Update tests/test_environment.py Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- tests/test_environment.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_environment.py b/tests/test_environment.py index db80eea6d..c0876248c 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -73,7 +73,7 @@ def test_wyoming_sounding_atmosphere(mock_show, example_env): pass assert example_env.all_info() == None assert abs(example_env.pressure(0) - 93600.0) < 1e-8 - assert example_env.barometric_height(example_env.pressure(0)) == 722.0 + assert abs(example_env.barometric_height(example_env.pressure(0)) - 722.0) < 1e-8 assert abs(example_env.wind_velocity_x(0) - -2.9005178894925043) < 1e-8 assert abs(example_env.temperature(100) - 291.75) < 1e-8 From 5349743e136eae1989c85697497b45891768c990 Mon Sep 17 00:00:00 2001 From: MateusStano <69485049+MateusStano@users.noreply.github.com> Date: Thu, 25 Jan 2024 12:36:29 -0300 Subject: [PATCH 20/25] Update tests/test_environment.py Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- tests/test_environment.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_environment.py b/tests/test_environment.py index c0876248c..608fe7e90 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -44,7 +44,7 @@ def test_custom_atmosphere(mock_show, example_env): ) assert example_env.all_info() == None assert abs(example_env.pressure(0) - 101325.0) < 1e-8 - assert example_env.barometric_height(101325.0) < 1e-2 + assert abs(example_env.barometric_height(101325.0)) < 1e-2 assert abs(example_env.wind_velocity_x(0) - 5) < 1e-8 assert abs(example_env.temperature(100) - 300) < 1e-8 From 09d0b4be2eaa423ff5b92178258624a21351ec8d Mon Sep 17 00:00:00 2001 From: MateusStano <69485049+MateusStano@users.noreply.github.com> Date: Thu, 25 Jan 2024 12:38:00 -0300 Subject: [PATCH 21/25] Update tests/test_environment.py Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- tests/test_environment.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_environment.py b/tests/test_environment.py index 608fe7e90..3067530b7 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -20,7 +20,7 @@ def test_standard_atmosphere(mock_show, example_env): assert example_env.info() == None assert example_env.all_info() == None assert abs(example_env.pressure(0) - 101325.0) < 1e-8 - assert example_env.barometric_height(101325.0) < 1e-2 + assert abs(example_env.barometric_height(101325.0)) < 1e-2 assert example_env.prints.print_earth_details() == None From 0de47560d35ae2cfe7da292fad3116a29b58e309 Mon Sep 17 00:00:00 2001 From: MateusStano <69485049+MateusStano@users.noreply.github.com> Date: Thu, 25 Jan 2024 12:38:22 -0300 Subject: [PATCH 22/25] Update rocketpy/environment/environment.py Co-authored-by: Gui-FernandesBR <63590233+Gui-FernandesBR@users.noreply.github.com> --- rocketpy/environment/environment.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 2778f19fb..8b8d0498a 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -3025,9 +3025,7 @@ def pressure_function(h): outputs="Pressure (Pa)", ) - # Save international standard atmosphere height by pressure profile - # and discretize it. This is done to speed up the calculations in the - # trajectory simulation. + # Discretize Function to speed up the trajectory simulation. self.barometric_height_ISA = self.pressure_ISA.inverse_function().set_discrete( pressure[-1], pressure[0], 100, extrapolation="constant" ) From b469931880dc574c81e7e252e55160c7cf5ad91a Mon Sep 17 00:00:00 2001 From: Lucas de Oliveira Prates Date: Thu, 25 Jan 2024 14:09:20 -0300 Subject: [PATCH 23/25] doc: Including imageio in pip install in dispersion analysis notebook. --- docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb b/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb index a15dc210f..a7794f7d5 100644 --- a/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb +++ b/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb @@ -66,7 +66,7 @@ }, "outputs": [], "source": [ - "!pip install rocketpy" + "!pip install rocketpy imageio" ] }, { From 0d57206ad42de34d8bcab100611a7bdca7d4f6a9 Mon Sep 17 00:00:00 2001 From: Lucas de Oliveira Prates Date: Wed, 31 Jan 2024 12:08:39 -0300 Subject: [PATCH 24/25] DOC: Replacing git clone command with curl in notebooks. --- .../deployable_payload_example.ipynb | 30 +++++++------------ .../dispersion_analysis.ipynb | 14 +++++---- .../parachute_drop_from_helicopter.ipynb | 14 +++++---- docs/notebooks/getting_started_colab.ipynb | 30 ++++++------------- 4 files changed, 36 insertions(+), 52 deletions(-) diff --git a/docs/notebooks/deployable_payload_example.ipynb b/docs/notebooks/deployable_payload_example.ipynb index fb96eb06d..38549ea72 100644 --- a/docs/notebooks/deployable_payload_example.ipynb +++ b/docs/notebooks/deployable_payload_example.ipynb @@ -33,7 +33,7 @@ "\n", "* RocketPy\n", "* netCDF4 (to get weather forecasts)\n", - "* Data files (we will clone RocketPy's repository for these)\n", + "* Data files (we will download RocketPy's data from github)\n", "\n", "Therefore, let's run the following lines of code:" ] @@ -45,18 +45,10 @@ "outputs": [], "source": [ "!pip install rocketpy netCDF4\n", - "!git clone https://github.com/RocketPy-Team/RocketPy.git" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "\n", - "os.chdir(\"RocketPy/docs/notebooks\")" + "!curl -o NACA0012-radians.csv https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/calisto/NACA0012-radians.csv\n", + "!curl -o Cesaroni_M1670.eng https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/motors/Cesaroni_M1670.eng\n", + "!curl -o powerOffDragCurve.csv https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/calisto/powerOffDragCurve.csv\n", + "!curl -o powerOnDragCurve.csv https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/calisto/powerOnDragCurve.csv" ] }, { @@ -237,7 +229,7 @@ "outputs": [], "source": [ "Pro75M1670 = SolidMotor(\n", - " thrust_source=\"../../data/motors/Cesaroni_M1670.eng\",\n", + " thrust_source=\"Cesaroni_M1670.eng\",\n", " dry_mass=1.815,\n", " dry_inertia=(0.125, 0.125, 0.002),\n", " nozzle_radius=33 / 1000,\n", @@ -392,8 +384,8 @@ " radius=127 / 2000,\n", " mass=rocket_mass + rocket_mass,\n", " inertia=(6.321, 6.321, 0.034),\n", - " power_off_drag=\"../../data/calisto/powerOffDragCurve.csv\",\n", - " power_on_drag=\"../../data/calisto/powerOnDragCurve.csv\",\n", + " power_off_drag=\"powerOffDragCurve.csv\",\n", + " power_on_drag=\"powerOnDragCurve.csv\",\n", " center_of_mass_without_motor=0,\n", " coordinate_system_orientation=\"tail_to_nose\",\n", ")\n", @@ -417,7 +409,7 @@ " span=0.110,\n", " position=-1.04956,\n", " cant_angle=0.5,\n", - " airfoil=(\"../../data/calisto/NACA0012-radians.csv\",\"radians\"),\n", + " airfoil=(\"NACA0012-radians.csv\",\"radians\"),\n", ")\n", "\n", "rocket_with_payload.add_tail(\n", @@ -536,8 +528,8 @@ " radius=127 / 2000,\n", " mass=rocket_mass,\n", " inertia=(6.321, 6.321, 0.034),\n", - " power_off_drag=\"../../data/calisto/powerOffDragCurve.csv\",\n", - " power_on_drag=\"../../data/calisto/powerOnDragCurve.csv\",\n", + " power_off_drag=\"powerOffDragCurve.csv\",\n", + " power_on_drag=\"powerOnDragCurve.csv\",\n", " center_of_mass_without_motor=0,\n", " coordinate_system_orientation=\"tail_to_nose\",\n", ")\n", diff --git a/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb b/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb index a7794f7d5..418589a28 100644 --- a/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb +++ b/docs/notebooks/dispersion_analysis/dispersion_analysis.ipynb @@ -25,10 +25,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Clone repository if using Google Colab\n", + "## Download data files if using Google Colab\n", "\n", "If you are running this using Binder, or you are running locally with the necessary files, you do not need to run this.\n", - "On the other hand, if you are running on Google Colab, make sure to run the cell below to clone the repository and download the necessary files." + "On the other hand, if you are running on Google Colab, make sure to run the cell below to download the necessary files." ] }, { @@ -37,10 +37,12 @@ "metadata": {}, "outputs": [], "source": [ - "!git clone https://github.com/RocketPy-Team/RocketPy.git\n", - "import os\n", - "\n", - "os.chdir(\"RocketPy/docs/notebooks/dispersion_analysis\")" + "!curl -o dispersion_analysis_inputs/Cd_PowerOff.csv --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/Cd_PowerOff.csv\n", + "!curl -o dispersion_analysis_inputs/Cd_PowerOn.csv --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/Cd_PowerOn.csv\n", + "!curl -o dispersion_analysis_inputs/LASC2019_reanalysis.nc --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/LASC2019_reanalysis.nc\n", + "!curl -o dispersion_analysis_inputs/thrustCurve.csv --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/thrustCurve.csv\n", + "!curl -o dispersion_analysis_inputs/Valetudo_basemap_final.jpg --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/Valetudo_basemap_final.jpg\n", + "!mkdir -p dispersion_analysis_outputs" ] }, { diff --git a/docs/notebooks/dispersion_analysis/parachute_drop_from_helicopter.ipynb b/docs/notebooks/dispersion_analysis/parachute_drop_from_helicopter.ipynb index e9709ff14..fd9e34663 100644 --- a/docs/notebooks/dispersion_analysis/parachute_drop_from_helicopter.ipynb +++ b/docs/notebooks/dispersion_analysis/parachute_drop_from_helicopter.ipynb @@ -19,10 +19,10 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Clone repository if using Google Colab\n", + "## Clone data files if using Google Colab\n", "\n", "If you are running this using Binder, or you are running locally with the necessary files, you do not need to run this.\n", - "On the other hand, if you are running on Google Colab, make sure to run the cell below to clone the repository and download the necessary files." + "On the other hand, if you are running on Google Colab, make sure to run the cell below to download the necessary files." ] }, { @@ -31,10 +31,12 @@ "metadata": {}, "outputs": [], "source": [ - "!git clone https://github.com/giovaniceotto/RocketPy.git\n", - "import os\n", - "\n", - "os.chdir(\"RocketPy/docs/notebooks/dispersion_analysis\")" + "!curl -o dispersion_analysis_inputs/Cd_PowerOff.csv --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/Cd_PowerOff.csv\n", + "!curl -o dispersion_analysis_inputs/Cd_PowerOn.csv --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/Cd_PowerOn.csv\n", + "!curl -o dispersion_analysis_inputs/LASC2019_reanalysis.nc --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/LASC2019_reanalysis.nc\n", + "!curl -o dispersion_analysis_inputs/thrustCurve.csv --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/thrustCurve.csv\n", + "!curl -o dispersion_analysis_inputs/Valetudo_basemap_final.jpg --create-dirs https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/docs/notebooks/dispersion_analysis/dispersion_analysis_inputs/Valetudo_basemap_final.jpg\n", + "!mkdir -p dispersion_analysis_outputs" ] }, { diff --git a/docs/notebooks/getting_started_colab.ipynb b/docs/notebooks/getting_started_colab.ipynb index 4ba05dd70..54f267ad7 100644 --- a/docs/notebooks/getting_started_colab.ipynb +++ b/docs/notebooks/getting_started_colab.ipynb @@ -22,7 +22,7 @@ "We start by setting up our environment. To run this notebook, we will need:\n", "\n", "- RocketPy\n", - "- Data files (we will clone RocketPy's repository for these)\n", + "- Data files (we will download RocketPy's data from github)\n", "\n", "Therefore, let's run the following lines of code:" ] @@ -38,22 +38,10 @@ "outputs": [], "source": [ "!pip install rocketpy\n", - "!git clone https://github.com/giovaniceotto/RocketPy.git" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "colab": {}, - "colab_type": "code", - "id": "pY5XGge5OoGJ" - }, - "outputs": [], - "source": [ - "import os\n", - "\n", - "os.chdir(\"RocketPy/docs/notebooks\")" + "!curl -o NACA0012-radians.csv https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/calisto/NACA0012-radians.csv\n", + "!curl -o Cesaroni_M1670.eng https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/motors/Cesaroni_M1670.eng\n", + "!curl -o powerOffDragCurve.csv https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/calisto/powerOffDragCurve.csv\n", + "!curl -o powerOnDragCurve.csv https://raw.githubusercontent.com/RocketPy-Team/RocketPy/master/data/calisto/powerOnDragCurve.csv" ] }, { @@ -2485,7 +2473,7 @@ "outputs": [], "source": [ "Pro75M1670 = SolidMotor(\n", - " thrust_source=\"../../data/motors/Cesaroni_M1670.eng\",\n", + " thrust_source=\"Cesaroni_M1670.eng\",\n", " dry_mass=1.815,\n", " dry_inertia=(0.125, 0.125, 0.002),\n", " nozzle_radius=33 / 1000,\n", @@ -3499,8 +3487,8 @@ " radius=127 / 2000,\n", " mass=14.426,\n", " inertia=(6.321, 6.321, 0.034),\n", - " power_off_drag=\"../../data/calisto/powerOffDragCurve.csv\",\n", - " power_on_drag=\"../../data/calisto/powerOnDragCurve.csv\",\n", + " power_off_drag=\"powerOffDragCurve.csv\",\n", + " power_on_drag=\"powerOnDragCurve.csv\",\n", " center_of_mass_without_motor=0,\n", " coordinate_system_orientation=\"tail_to_nose\",\n", ")\n", @@ -3579,7 +3567,7 @@ " span=0.110,\n", " position=-1.04956,\n", " cant_angle=0.5,\n", - " airfoil=(\"../../data/calisto/NACA0012-radians.csv\",\"radians\"),\n", + " airfoil=(\"NACA0012-radians.csv\",\"radians\"),\n", ")\n", "\n", "tail = calisto.add_tail(\n", From be46ce8b6d56697e0fcfe05b342e5d92c231e685 Mon Sep 17 00:00:00 2001 From: Lint Action Date: Wed, 31 Jan 2024 15:38:37 +0000 Subject: [PATCH 25/25] Fix code style issues with Black --- rocketpy/environment/environment_analysis.py | 8 +-- rocketpy/mathutils/function.py | 1 + rocketpy/plots/flight_plots.py | 40 ++++++++----- rocketpy/rocket/aero_surface.py | 8 +-- rocketpy/simulation/flight.py | 59 +++++++++----------- rocketpy/simulation/flight_data_importer.py | 1 + tests/test_flight_data_importer.py | 1 + tests/test_genericmotor.py | 4 +- tests/test_tank.py | 11 +--- tests/unit/test_rocket.py | 3 +- tests/unit/test_solidmotor.py | 4 +- 11 files changed, 65 insertions(+), 75 deletions(-) diff --git a/rocketpy/environment/environment_analysis.py b/rocketpy/environment/environment_analysis.py index 2be707d73..c15b32551 100644 --- a/rocketpy/environment/environment_analysis.py +++ b/rocketpy/environment/environment_analysis.py @@ -901,10 +901,10 @@ def __parse_surface_data(self): # Extract data from weather file indices = (time_index, lon_index, lat_index) for key, value in surface_file_dict.items(): - dictionary[date_string][hour_string][ - key - ] = self.__extract_surface_data_value( - surface_data, value, indices, lon_array, lat_array + dictionary[date_string][hour_string][key] = ( + self.__extract_surface_data_value( + surface_data, value, indices, lon_array, lat_array + ) ) # Get elevation, time index does not matter, use last one diff --git a/rocketpy/mathutils/function.py b/rocketpy/mathutils/function.py index 7981bbbfa..c899f3cad 100644 --- a/rocketpy/mathutils/function.py +++ b/rocketpy/mathutils/function.py @@ -3,6 +3,7 @@ and more. This is a core class of our package, and should be maintained carefully as it may impact all the rest of the project. """ + import warnings from collections.abc import Iterable from copy import deepcopy diff --git a/rocketpy/plots/flight_plots.py b/rocketpy/plots/flight_plots.py index e3b144d97..289daf5eb 100644 --- a/rocketpy/plots/flight_plots.py +++ b/rocketpy/plots/flight_plots.py @@ -400,9 +400,11 @@ def rail_buttons_forces(self): ) ax1.set_xlim( 0, - self.flight.out_of_rail_time - if self.flight.out_of_rail_time > 0 - else self.flight.tFinal, + ( + self.flight.out_of_rail_time + if self.flight.out_of_rail_time > 0 + else self.flight.tFinal + ), ) ax1.legend() ax1.grid(True) @@ -431,9 +433,11 @@ def rail_buttons_forces(self): ) ax2.set_xlim( 0, - self.flight.out_of_rail_time - if self.flight.out_of_rail_time > 0 - else self.flight.tFinal, + ( + self.flight.out_of_rail_time + if self.flight.out_of_rail_time > 0 + else self.flight.tFinal + ), ) ax2.legend() ax2.grid(True) @@ -557,9 +561,11 @@ def energy_data(self): ) ax1.set_xlim( 0, - self.flight.apogee_time - if self.flight.apogee_time != 0.0 - else self.flight.t_final, + ( + self.flight.apogee_time + if self.flight.apogee_time != 0.0 + else self.flight.t_final + ), ) ax1.ticklabel_format(style="sci", axis="y", scilimits=(0, 0)) ax1.set_title("Kinetic Energy Components") @@ -587,9 +593,11 @@ def energy_data(self): ) ax2.set_xlim( 0, - self.flight.apogee_time - if self.flight.apogee_time != 0.0 - else self.flight.t_final, + ( + self.flight.apogee_time + if self.flight.apogee_time != 0.0 + else self.flight.t_final + ), ) ax2.ticklabel_format(style="sci", axis="y", scilimits=(0, 0)) ax2.set_title("Total Mechanical Energy Components") @@ -620,9 +628,11 @@ def energy_data(self): ) ax4.set_xlim( 0, - self.flight.apogee_time - if self.flight.apogee_time != 0.0 - else self.flight.t_final, + ( + self.flight.apogee_time + if self.flight.apogee_time != 0.0 + else self.flight.t_final + ), ) ax3.ticklabel_format(style="sci", axis="y", scilimits=(0, 0)) ax4.set_title("Drag Absolute Power") diff --git a/rocketpy/rocket/aero_surface.py b/rocketpy/rocket/aero_surface.py index cea919141..a72fb74cc 100644 --- a/rocketpy/rocket/aero_surface.py +++ b/rocketpy/rocket/aero_surface.py @@ -1180,10 +1180,7 @@ def evaluate_geometrical_parameters(self): * (self.root_chord + 2 * self.tip_chord) * self.rocket_radius * self.span**2 - + 6 - * (self.root_chord + self.tip_chord) - * self.span - * self.rocket_radius**2 + + 6 * (self.root_chord + self.tip_chord) * self.span * self.rocket_radius**2 ) / 12 roll_damping_interference_factor = 1 + ( ((tau - λ) / (tau)) - ((1 - λ) / (tau - 1)) * np.log(tau) @@ -1506,8 +1503,7 @@ def evaluate_geometrical_parameters(self): * self.rocket_radius**2 * np.sqrt(-self.span**2 + self.rocket_radius**2) * np.arctan( - (self.span) - / (np.sqrt(-self.span**2 + self.rocket_radius**2)) + (self.span) / (np.sqrt(-self.span**2 + self.rocket_radius**2)) ) - np.pi * self.rocket_radius**2 diff --git a/rocketpy/simulation/flight.py b/rocketpy/simulation/flight.py index 39b0c5c7e..d72370140 100644 --- a/rocketpy/simulation/flight.py +++ b/rocketpy/simulation/flight.py @@ -683,9 +683,9 @@ def __init__( node.time_bound = phase.TimeNodes[node_index + 1].t phase.solver.t_bound = node.time_bound phase.solver._lsoda_solver._integrator.rwork[0] = phase.solver.t_bound - phase.solver._lsoda_solver._integrator.call_args[ - 4 - ] = phase.solver._lsoda_solver._integrator.rwork + phase.solver._lsoda_solver._integrator.call_args[4] = ( + phase.solver._lsoda_solver._integrator.rwork + ) phase.solver.status = "running" # Feed required parachute and discrete controller triggers @@ -1275,9 +1275,7 @@ def udot_rail1(self, t, u, post_processing=False): R3 = -0.5 * rho * (free_stream_speed**2) * self.rocket.area * (drag_coeff) # Calculate Linear acceleration - a3 = (R3 + thrust) / M - ( - e0**2 - e1**2 - e2**2 + e3**2 - ) * self.env.gravity(z) + a3 = (R3 + thrust) / M - (e0**2 - e1**2 - e2**2 + e3**2) * self.env.gravity(z) if a3 > 0: ax = 2 * (e1 * e3 + e0 * e2) * a3 ay = 2 * (e2 * e3 - e0 * e1) * a3 @@ -1469,9 +1467,7 @@ def u_dot(self, t, u, post_processing=False): 0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift ) # component lift force components - lift_dir_norm = ( - comp_stream_vx_b**2 + comp_stream_vy_b**2 - ) ** 0.5 + lift_dir_norm = (comp_stream_vx_b**2 + comp_stream_vy_b**2) ** 0.5 comp_lift_xb = comp_lift * (comp_stream_vx_b / lift_dir_norm) comp_lift_yb = comp_lift * (comp_stream_vy_b / lift_dir_norm) # add to total lift force @@ -1750,9 +1746,7 @@ def u_dot_generalized(self, t, u, post_processing=False): 0.5 * rho * (comp_stream_speed**2) * reference_area * c_lift ) # Component lift force components - lift_dir_norm = ( - comp_stream_vx_b**2 + comp_stream_vy_b**2 - ) ** 0.5 + lift_dir_norm = (comp_stream_vx_b**2 + comp_stream_vy_b**2) ** 0.5 comp_lift_xb = comp_lift * (comp_stream_vx_b / lift_dir_norm) comp_lift_yb = comp_lift * (comp_stream_vy_b / lift_dir_norm) # Add to total lift force @@ -1896,8 +1890,7 @@ def u_dot_parachute(self, t, u, post_processing=False): freestream_z = vz # Determine drag force pseudoD = ( - -0.5 * rho * cd_s * free_stream_speed - - ka * rho * 4 * np.pi * (R**2) * Rdot + -0.5 * rho * cd_s * free_stream_speed - ka * rho * 4 * np.pi * (R**2) * Rdot ) Dx = pseudoD * freestream_x Dy = pseudoD * freestream_y @@ -2486,9 +2479,7 @@ def translational_energy(self): # Redefine total_mass time grid to allow for efficient Function algebra total_mass = deepcopy(self.rocket.total_mass) total_mass.set_discrete_based_on_model(self.vz) - translational_energy = ( - 0.5 * total_mass * (self.vx**2 + self.vy**2 + self.vz**2) - ) + translational_energy = 0.5 * total_mass * (self.vx**2 + self.vy**2 + self.vz**2) return translational_energy @funcify_method("Time (s)", "Kinetic Energy (J)", "spline", "zero") @@ -3400,19 +3391,23 @@ def add(self, flight_phase, index=None): ) if flight_phase.t == previous_phase.t else ( - "Trying to add flight phase starting *together* with the one *proceeding* it. ", - "This may be caused by multiple parachutes being triggered simultaneously.", - ) - if flight_phase.t == next_phase.t - else ( - "Trying to add flight phase starting *before* the one *preceding* it. ", - "This may be caused by multiple parachutes being triggered simultaneously", - " or by having a negative parachute lag.", - ) - if flight_phase.t < previous_phase.t - else ( - "Trying to add flight phase starting *after* the one *proceeding* it.", - "This may be caused by multiple parachutes being triggered simultaneously.", + ( + "Trying to add flight phase starting *together* with the one *proceeding* it. ", + "This may be caused by multiple parachutes being triggered simultaneously.", + ) + if flight_phase.t == next_phase.t + else ( + ( + "Trying to add flight phase starting *before* the one *preceding* it. ", + "This may be caused by multiple parachutes being triggered simultaneously", + " or by having a negative parachute lag.", + ) + if flight_phase.t < previous_phase.t + else ( + "Trying to add flight phase starting *after* the one *proceeding* it.", + "This may be caused by multiple parachutes being triggered simultaneously.", + ) + ) ) ) self.display_warning(*warning_msg) @@ -3420,9 +3415,7 @@ def add(self, flight_phase, index=None): new_index = ( index - 1 if flight_phase.t < previous_phase.t - else index + 1 - if flight_phase.t > next_phase.t - else index + else index + 1 if flight_phase.t > next_phase.t else index ) flight_phase.t += adjust self.add(flight_phase, new_index) diff --git a/rocketpy/simulation/flight_data_importer.py b/rocketpy/simulation/flight_data_importer.py index 3b6d3f473..edb2390d6 100644 --- a/rocketpy/simulation/flight_data_importer.py +++ b/rocketpy/simulation/flight_data_importer.py @@ -1,6 +1,7 @@ """Starts with .csv or .txt file containing the rocket's collected flight data and build a rocketpy.Flight object from it. """ + import warnings from os import listdir from os.path import isfile, join diff --git a/tests/test_flight_data_importer.py b/tests/test_flight_data_importer.py index 9ce4e34ce..2dd66f3ae 100644 --- a/tests/test_flight_data_importer.py +++ b/tests/test_flight_data_importer.py @@ -1,5 +1,6 @@ """Tests the FlightDataImporter class from rocketpy.simulation module. """ + import numpy as np from rocketpy.simulation import FlightDataImporter diff --git a/tests/test_genericmotor.py b/tests/test_genericmotor.py index d9c3ff5fd..513fca40d 100644 --- a/tests/test_genericmotor.py +++ b/tests/test_genericmotor.py @@ -118,9 +118,7 @@ def test_generic_motor_inertia(generic_motor): # Tests the inertia formulation from the propellant mass propellant_mass = generic_motor.propellant_mass.set_discrete(*burn_time, 50).y_array - propellant_I_11 = propellant_mass * ( - chamber_radius**2 / 4 + chamber_height**2 / 12 - ) + propellant_I_11 = propellant_mass * (chamber_radius**2 / 4 + chamber_height**2 / 12) propellant_I_22 = propellant_I_11 propellant_I_33 = propellant_mass * (chamber_radius**2 / 2) diff --git a/tests/test_tank.py b/tests/test_tank.py index 9a6ddb6b0..6b2f03bc9 100644 --- a/tests/test_tank.py +++ b/tests/test_tank.py @@ -631,11 +631,7 @@ def upper_spherical_cap_centroid(radius, height=None): """ if height is None: height = radius - return ( - 0.75 - * (height**3 - 2 * height * radius**2) - / (height**2 - 3 * radius**2) - ) + return 0.75 * (height**3 - 2 * height * radius**2) / (height**2 - 3 * radius**2) def tank_centroid_function(tank_radius, tank_height, zero_height=0): @@ -775,10 +771,7 @@ def lower_spherical_cap_inertia(radius, height=None, reference=0): ) inertia_y = inertia_x inertia_z = lower_spherical_cap_volume(radius, height) * ( - np.pi - * height**3 - * (3 * height**2 - 15 * height * radius + 20 * radius**2) - / 30 + np.pi * height**3 * (3 * height**2 - 15 * height * radius + 20 * radius**2) / 30 ) return np.array([inertia_x, inertia_y, inertia_z]) diff --git a/tests/unit/test_rocket.py b/tests/unit/test_rocket.py index 1d8cfac21..7d942ee60 100644 --- a/tests/unit/test_rocket.py +++ b/tests/unit/test_rocket.py @@ -207,8 +207,7 @@ def test_add_fins_assert_cp_cm_plus_fins(calisto, dimensionless_calisto, m): 1 + np.sqrt( 1 - + (2 * np.sqrt((0.12 / 2 - 0.04 / 2) ** 2 + 0.1**2) / (0.120 + 0.040)) - ** 2 + + (2 * np.sqrt((0.12 / 2 - 0.04 / 2) ** 2 + 0.1**2) / (0.120 + 0.040)) ** 2 ) ) clalpha *= 1 + calisto.radius / (0.1 + calisto.radius) diff --git a/tests/unit/test_solidmotor.py b/tests/unit/test_solidmotor.py index b23dfc31e..6ae0f68ca 100644 --- a/tests/unit/test_solidmotor.py +++ b/tests/unit/test_solidmotor.py @@ -70,9 +70,7 @@ def test_evaluate_inertia_33_asserts_extreme_values(cesaroni_m1670): grain_mass = grain_vol * grain_density grain_I_33_initial = ( - grain_mass - * (1 / 2.0) - * (grain_initial_inner_radius**2 + grain_outer_radius**2) + grain_mass * (1 / 2.0) * (grain_initial_inner_radius**2 + grain_outer_radius**2) ) # not passing because I_33 is not discrete anymore