From 73097568139be65fabd5b24c6c57a6e95f622f00 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sat, 16 Sep 2023 14:52:38 -0300 Subject: [PATCH 1/4] DEP: remove deprecated or useless methods --- rocketpy/Flight.py | 63 ---------- rocketpy/Function.py | 255 ----------------------------------------- tests/test_function.py | 2 - 3 files changed, 320 deletions(-) diff --git a/rocketpy/Flight.py b/rocketpy/Flight.py index f1e9fa07a..91f49b555 100644 --- a/rocketpy/Flight.py +++ b/rocketpy/Flight.py @@ -1,10 +1,8 @@ import math -import time import warnings from copy import deepcopy from functools import cached_property -import matplotlib.pyplot as plt import numpy as np import simplekml from scipy import integrate @@ -3313,67 +3311,6 @@ def all_info(self): return None - def animate(self, start=0, stop=None, fps=12, speed=4, elev=None, azim=None): - """Plays an animation of the flight. Not implemented yet. Only - kinda works outside notebook. - """ - # Set up stopping time - stop = self.t_final if stop is None else stop - # Speed = 4 makes it almost real time - matplotlib is way to slow - # Set up graph - fig = plt.figure(figsize=(18, 15)) - axes = fig.gca(projection="3d") - # Initialize time - time_range = np.linspace(start, stop, fps * (stop - start)) - # Initialize first frame - axes.set_title("Trajectory and Velocity Animation") - axes.set_xlabel("X (m)") - axes.set_ylabel("Y (m)") - axes.set_zlabel("Z (m)") - axes.view_init(elev, azim) - R = axes.quiver(0, 0, 0, 0, 0, 0, color="r", label="Rocket") - V = axes.quiver(0, 0, 0, 0, 0, 0, color="g", label="Velocity") - W = axes.quiver(0, 0, 0, 0, 0, 0, color="b", label="Wind") - S = axes.quiver(0, 0, 0, 0, 0, 0, color="black", label="Freestream") - axes.legend() - # Animate - for t in time_range: - R.remove() - V.remove() - W.remove() - S.remove() - # Calculate rocket position - Rx, Ry, Rz = self.x(t), self.y(t), self.z(t) - Ru = 1 * (2 * (self.e1(t) * self.e3(t) + self.e0(t) * self.e2(t))) - Rv = 1 * (2 * (self.e2(t) * self.e3(t) - self.e0(t) * self.e1(t))) - Rw = 1 * (1 - 2 * (self.e1(t) ** 2 + self.e2(t) ** 2)) - # Calculate rocket Mach number - Vx = self.vx(t) / 340.40 - Vy = self.vy(t) / 340.40 - Vz = self.vz(t) / 340.40 - # Calculate wind Mach Number - z = self.z(t) - Wx = self.env.wind_velocity_x(z) / 20 - Wy = self.env.wind_velocity_y(z) / 20 - # Calculate freestream Mach Number - Sx = self.stream_velocity_x(t) / 340.40 - Sy = self.stream_velocity_y(t) / 340.40 - Sz = self.stream_velocity_z(t) / 340.40 - # Plot Quivers - R = axes.quiver(Rx, Ry, Rz, Ru, Rv, Rw, color="r") - V = axes.quiver(Rx, Ry, Rz, -Vx, -Vy, -Vz, color="g") - W = axes.quiver(Rx - Vx, Ry - Vy, Rz - Vz, Wx, Wy, 0, color="b") - S = axes.quiver(Rx, Ry, Rz, Sx, Sy, Sz, color="black") - # Adjust axis - axes.set_xlim(Rx - 1, Rx + 1) - axes.set_ylim(Ry - 1, Ry + 1) - axes.set_zlim(Rz - 1, Rz + 1) - # plt.pause(1/(fps*speed)) - try: - plt.pause(1 / (fps * speed)) - except: - time.sleep(1 / (fps * speed)) - def time_iterator(self, node_list): i = 0 while i < len(node_list) - 1: diff --git a/rocketpy/Function.py b/rocketpy/Function.py index 291cf4fd9..7209ef4be 100644 --- a/rocketpy/Function.py +++ b/rocketpy/Function.py @@ -885,261 +885,6 @@ def get_value(self, *args): else: return x if len(x) > 1 else x[0] - def get_value_opt_deprecated(self, *args): - """THE CODE BELOW IS HERE FOR DOCUMENTATION PURPOSES ONLY. IT WAS - REPLACED FOR ALL INSTANCES BY THE FUNCTION.SETGETVALUEOPT METHOD. - - This method returns the value of the Function at the specified - point in a limited but optimized manner. See Function.get_value for an - implementation which allows more kinds of inputs. - This method optimizes the Function.get_value method by only - implementing function evaluations of single inputs, i.e., it is not - vectorized. Furthermore, it actually implements a different method - for each interpolation type, eliminating some if statements. - Currently supports callables and spline, linear, akima, polynomial and - shepard interpolated Function objects. - - Parameters - ---------- - args : scalar - Value where the Function is to be evaluated. If the Function is - 1-D, only one argument is expected, which may be an int or a float - If the function is N-D, N arguments must be given, each one being - an int or a float. - - Returns - ------- - x : scalar - """ - # Callables - if callable(self.source): - return self.source(*args) - - # Interpolated Function - # Retrieve general info - x_data = self.x_array - y_data = self.y_array - x_min, x_max = self.xinitial, self.xfinal - if self.__extrapolation__ == "zero": - extrapolation = 0 # Extrapolation is zero - elif self.__extrapolation__ == "natural": - extrapolation = 1 # Extrapolation is natural - elif self.__extrapolation__ == "constant": - extrapolation = 2 # Extrapolation is constant - else: - raise ValueError(f"Invalid extrapolation type {self.__extrapolation__}") - - # Interpolate this info for each interpolation type - # Spline - if self.__interpolation__ == "spline": - x = args[0] - coeffs = self.__spline_coefficients__ - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - x_interval = x_interval if x_interval != 0 else 1 - a = coeffs[:, x_interval - 1] - x = x - x_data[x_interval - 1] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - a = coeffs[:, 0] if x < x_min else coeffs[:, -1] - x = x - x_data[0] if x < x_min else x - x_data[-2] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - # Linear - elif self.__interpolation__ == "linear": - x = args[0] - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - y = (x - x_data[x_interval - 1]) * (dy / dx) + y_data[x_interval - 1] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - x_interval = 1 if x < x_min else -1 - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - y = (x - x_data[x_interval - 1]) * (dy / dx) + y_data[ - x_interval - 1 - ] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - # Akima - elif self.__interpolation__ == "akima": - x = args[0] - coeffs = np.array(self.__akima_coefficients__) - x_interval = np.searchsorted(x_data, x) - # Interval found... interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - x_interval = x_interval if x_interval != 0 else 1 - a = coeffs[4 * x_interval - 4 : 4 * x_interval] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - a = coeffs[:4] if x < x_min else coeffs[-4:] - y = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - # Polynomial - elif self.__interpolation__ == "polynomial": - x = args[0] - coeffs = self.__polynomial_coefficients__ - # Interpolate... or extrapolate - if x_min <= x <= x_max: - # Interpolate - y = 0 - for i in range(len(coeffs)): - y += coeffs[i] * (x**i) - else: - # Extrapolate - if extrapolation == 0: # Extrapolation == zero - y = 0 - elif extrapolation == 1: # Extrapolation == natural - y = 0 - for i in range(len(coeffs)): - y += coeffs[i] * (x**i) - else: # Extrapolation is set to constant - y = y_data[0] if x < x_min else y_data[-1] - return y - # Shepard - elif self.__interpolation__ == "shepard": - x_data = self.source[:, 0:-1] # Support for N-Dimensions - len_y_data = len(y_data) # A little speed up - x = np.array([[float(x) for x in list(args)]]) - numerator_sum = 0 - denominator_sum = 0 - for i in range(len_y_data): - sub = x_data[i] - x - distance = np.linalg.norm(sub) - if distance == 0: - numerator_sum = y_data[i] - denominator_sum = 1 - break - else: - weight = distance ** (-3) - numerator_sum = numerator_sum + y_data[i] * weight - denominator_sum = denominator_sum + weight - return numerator_sum / denominator_sum - - def get_value_opt2(self, *args): - """DEPRECATED!! - See Function.get_value_opt for new version. - This method returns the value of the Function at the specified - point in a limited but optimized manner. See Function.get_value for an - implementation which allows more kinds of inputs. - This method optimizes the Function.get_value method by only - implementing function evaluations of single inputs, i.e., it is not - vectorized. Furthermore, it actually implements a different method - for each interpolation type, eliminating some if statements. - Finally, it uses Numba to compile the methods, which further optimizes - the implementation. - The code below is here for documentation purposes only. It is - overwritten for all instances by the Function.setGetValuteOpt2 method. - - Parameters - ---------- - args : scalar - Value where the Function is to be evaluated. If the Function is - 1-D, only one argument is expected, which may be an int or a float - If the function is N-D, N arguments must be given, each one being - an int or a float. - - Returns - ------- - x : scalar - """ - # Returns value for function function type - if callable(self.source): - return self.source(*args) - # Returns value for spline, akima or linear interpolation function type - elif self.__interpolation__ in ["spline", "akima", "linear"]: - x = args[0] - x_data = self.x_array - y_data = self.y_array - # Hunt in intervals near the last interval which was used. - x_interval = self.last_interval - if x_data[x_interval - 1] <= x <= x_data[x_interval]: - pass - else: - x_interval = np.searchsorted(x_data, x) - self.last_interval = x_interval if x_interval < len(x_data) else 0 - # Interval found... keep going - x_min, x_max = self.xinitial, self.xfinal - if self.__interpolation__ == "spline": - coeffs = self.__spline_coefficients__ - if x == x_min or x == x_max: - x = y_data[x_interval] - elif x_min < x < x_max or (self.__extrapolation__ == "natural"): - if not x_min < x < x_max: - a = coeffs[:, 0] if x < x_min else coeffs[:, -1] - x = x - x_data[0] if x < x_min else x - x_data[-2] - else: - a = coeffs[:, x_interval - 1] - x = x - x_data[x_interval - 1] - x = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if self.__extrapolation__ == "zero": - x = 0 - else: # Extrapolation is set to constant - x = y_data[0] if x < x_min else y_data[-1] - elif self.__interpolation__ == "linear": - if x == x_min or x == x_max: - x = y_data[x_interval] - elif x_min < x < x_max or (self.__extrapolation__ == "natural"): - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - x = (x - x_data[x_interval - 1]) * (dy / dx) + y_data[ - x_interval - 1 - ] - elif self.__extrapolation__ == "natural": - y0 = y_data[0] if x < x_min else y_data[-1] - x_interval = 1 if x < x_min else -1 - dx = float(x_data[x_interval] - x_data[x_interval - 1]) - dy = float(y_data[x_interval] - y_data[x_interval - 1]) - x = (x - x_data[x_interval - 1]) * (dy / dx) + y0 - else: - # Extrapolate - if self.__extrapolation__ == "zero": - x = 0 - else: # Extrapolation is set to constant - x = y_data[0] if x < x_min else y_data[-1] - else: - if self.__interpolation__ == "akima": - coeffs = self.__akima_coefficients__ - if x == x_min or x == x_max: - x = y_data[x_interval] - elif x_min < x < x_max: - a = coeffs[4 * x_interval - 4 : 4 * x_interval] - x = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - elif self.__extrapolation__ == "natural": - a = coeffs[:4] if x < x_min else coeffs[-4:] - x = a[3] * x**3 + a[2] * x**2 + a[1] * x + a[0] - else: - # Extrapolate - if self.__extrapolation__ == "zero": - x = 0 - else: # Extrapolation is set to constant - x = y_data[0] if x < x_min else y_data[-1] - return x - def __getitem__(self, args): """Returns item of the Function source. If the source is not an array, an error will result. diff --git a/tests/test_function.py b/tests/test_function.py index f989f4152..4cbac9c33 100644 --- a/tests/test_function.py +++ b/tests/test_function.py @@ -39,9 +39,7 @@ def test_getters(func_from_csv): assert func_from_csv.get_interpolation_method() == "linear" assert func_from_csv.get_extrapolation_method() == "natural" assert np.isclose(func_from_csv.get_value(0), 0.0, atol=1e-6) - assert np.isclose(func_from_csv.get_value_opt_deprecated(0), 0.0, atol=1e-6) assert np.isclose(func_from_csv.get_value_opt(0), 0.0, atol=1e-6) - assert np.isclose(func_from_csv.get_value_opt2(0), 0.0, atol=1e-6) def test_setters(func_from_csv): From 39b77b97f9d432404a683f5c3d1ff184a2bad2a1 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sat, 16 Sep 2023 14:53:18 -0300 Subject: [PATCH 2/4] BUG: conversion of decimal degrees not working - the problem was happening when the angle was negative. - i.e. the operator // didn't work - I added some unit tests to cover the problematic cases Everything running properly now --- rocketpy/Environment.py | 29 +++++++++++------------------ tests/test_environment.py | 31 +++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 18 deletions(-) diff --git a/rocketpy/Environment.py b/rocketpy/Environment.py index e59380280..9912c1b06 100644 --- a/rocketpy/Environment.py +++ b/rocketpy/Environment.py @@ -3560,7 +3560,8 @@ def calculate_earth_radius(self, lat): return e_radius - def decimal_degrees_to_arc_seconds(self, angle): + @staticmethod + def decimal_degrees_to_arc_seconds(angle): """Function to convert an angle in decimal degrees to deg/min/sec. Converts (°) to (° ' ") @@ -3572,24 +3573,16 @@ def decimal_degrees_to_arc_seconds(self, angle): Returns ------- - deg : float + degrees : float The degrees. - min : float + arc_minutes : float The arc minutes. 1 arc-minute = (1/60)*degree - sec : float + arc_seconds : float The arc Seconds. 1 arc-second = (1/3600)*degree """ - - if angle < 0: - signal = -1 - else: - signal = 1 - - deg = (signal * angle) // 1 - min = abs(signal * angle - deg) * 60 // 1 - sec = abs((signal * angle - deg) * 60 - min) * 60 - # print("The angle {:f} is equals to {:.0f}º {:.0f}' {:.3f}'' ".format( - # angle, signal*deg, min, sec - # )) - - return deg, min, sec + sign = -1 if angle < 0 else 1 + degrees = int(abs(angle)) * sign + remainder = abs(angle) - abs(degrees) + arc_minutes = int(remainder * 60) + arc_seconds = (remainder * 60 - arc_minutes) * 60 + return degrees, arc_minutes, arc_seconds diff --git a/tests/test_environment.py b/tests/test_environment.py index 46d19a828..818e43303 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -7,6 +7,8 @@ import pytest import pytz +from rocketpy.Environment import Environment + def test_env_set_date(example_env): """Test that the date is set correctly in the environment object. This @@ -383,3 +385,32 @@ def test_utm_to_geodesic(example_env_robust): ) assert np.isclose(lat, 32.99025, atol=1e-5) == True assert np.isclose(lon, -106.9750, atol=1e-5) == True + +@pytest.mark.parametrize( + "angle, deg, arc_min, arc_sec", + [ + (-106.974998, -106.0, 58, 29.9928), + (32.990254, 32, 59.0, 24.9144), + (90.0, 90, 0, 0), + ], +) +def test_decimal_degrees_to_arc_seconds(angle, deg, arc_min, arc_sec): + """Tests if the conversion from decimal degrees to arc seconds is correct. + It takes 3 different angles and their expected results and compares them + with the results from the method. + + Parameters + ---------- + angle : float + Angle in decimal degrees. + deg : int + Expected degrees. + arc_min : int + Expected arc minutes. + arc_sec : float + Expected arc seconds. + """ + res = Environment.decimal_degrees_to_arc_seconds(angle) + assert pytest.approx(res[0], abs=1e-8) == deg + assert pytest.approx(res[1], abs=1e-8) == arc_min + assert pytest.approx(res[2], abs=1e-8) == arc_sec From d757c0b7fd65064982ca0d79bd86f855df0f9aa1 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Sat, 16 Sep 2023 14:53:59 -0300 Subject: [PATCH 3/4] MAINT: convert Environment methods to static --- rocketpy/Environment.py | 74 ++++++++++++++++++++++++++------------- tests/test_environment.py | 57 ++++++++++++++++++++++++------ 2 files changed, 97 insertions(+), 34 deletions(-) diff --git a/rocketpy/Environment.py b/rocketpy/Environment.py index 9912c1b06..913cd8a2f 100644 --- a/rocketpy/Environment.py +++ b/rocketpy/Environment.py @@ -362,7 +362,12 @@ def __init__( # Store launch site coordinates referenced to UTM projection system if self.latitude > -80 and self.latitude < 84: - convert = self.geodesic_to_utm(self.latitude, self.longitude) + convert = self.geodesic_to_utm( + lat=self.latitude, + lon=self.longitude, + flattening=self.ellipsoid.flattening, + semi_major_axis=self.ellipsoid.semi_major_axis, + ) self.initial_north = convert[1] self.initial_east = convert[0] @@ -375,8 +380,12 @@ def __init__( self.elevation = elevation self.set_elevation(elevation) - # Recalculate Earth Radius - self.earth_radius = self.calculate_earth_radius(self.latitude) # in m + # Recalculate Earth Radius (meters) + self.earth_radius = self.calculate_earth_radius( + lat=self.latitude, + semi_major_axis=self.ellipsoid.semi_major_axis, + flattening=self.ellipsoid.flattening, + ) # Initialize plots and prints object self.plots = _EnvironmentPlots(self) @@ -3309,7 +3318,11 @@ def set_earth_geometry(self, datum): ) # Auxiliary functions - Geodesic Coordinates - def geodesic_to_utm(self, lat, lon): + + @staticmethod + def geodesic_to_utm( + lat, lon, semi_major_axis=6378137.0, flattening=1 / 298.257223563 + ): """Function which converts geodetic coordinates, i.e. lat/lon, to UTM projection coordinates. Can be used only for latitudes between -80.00° and 84.00° @@ -3322,6 +3335,14 @@ def geodesic_to_utm(self, lat, lon): lon : float The longitude coordinates of the point of analysis, must be contained between -180.00° and 180.00° + semi_major_axis : float + The semi-major axis of the ellipsoid used to represent the Earth, + must be given in meters (default is 6,378,137.0 m, which corresponds + to the WGS84 ellipsoid) + flattening : float + The flattening of the ellipsoid used to represent the Earth, usually + between 1/250 and 1/150 (default is 1/298.257223563, which + corresponds to the WGS84 ellipsoid) Returns ------- @@ -3360,10 +3381,6 @@ def geodesic_to_utm(self, lat, lon): lon_mc = 3 EW = "W|E" - # Select the desired datum (i.e. the ellipsoid parameters) - flattening = self.ellipsoid.flattening - semi_major_axis = self.ellipsoid.semi_major_axis - # Evaluate the hemisphere and determine the N coordinate at the Equator if lat < 0: N0 = 10000000 @@ -3424,7 +3441,10 @@ def geodesic_to_utm(self, lat, lon): return x, y, utm_zone, utm_letter, hemis, EW - def utm_to_geodesic(self, x, y, utm_zone, hemis): + @staticmethod + def utm_to_geodesic( + x, y, utm_zone, hemis, semi_major_axis=6378137.0, flattening=1 / 298.257223563 + ): """Function to convert UTM coordinates to geodesic coordinates (i.e. latitude and longitude). The latitude should be between -80° and 84° @@ -3441,6 +3461,14 @@ def utm_to_geodesic(self, x, y, utm_zone, hemis): hemis : string Equals to "S" for southern hemisphere and "N" for Northern hemisphere + semi_major_axis : float + The semi-major axis of the ellipsoid used to represent the Earth, + must be given in meters (default is 6,378,137.0 m, which corresponds + to the WGS84 ellipsoid) + flattening : float + The flattening of the ellipsoid used to represent the Earth, usually + between 1/250 and 1/150 (default is 1/298.257223563, which + corresponds to the WGS84 ellipsoid) Returns ------- @@ -3456,10 +3484,6 @@ def utm_to_geodesic(self, x, y, utm_zone, hemis): # Calculate the Central Meridian from the UTM zone number central_meridian = utm_zone * 6 - 183 # degrees - # Select the desired datum - flattening = self.ellipsoid.flattening - semi_major_axis = self.ellipsoid.semi_major_axis - # Calculate reference values K0 = 1 - 1 / 2500 e2 = 2 * flattening - flattening**2 @@ -3513,7 +3537,10 @@ def utm_to_geodesic(self, x, y, utm_zone, hemis): return lat, lon - def calculate_earth_radius(self, lat): + @staticmethod + def calculate_earth_radius( + lat, semi_major_axis=6378137.0, flattening=1 / 298.257223563 + ): """Simple function to calculate the Earth Radius at a specific latitude based on ellipsoidal reference model (datum). The earth radius here is assumed as the distance between the ellipsoid's center of gravity and a @@ -3526,21 +3553,23 @@ def calculate_earth_radius(self, lat): ---------- lat : float latitude in which the Earth radius will be calculated + semi_major_axis : float + The semi-major axis of the ellipsoid used to represent the Earth, + must be given in meters (default is 6,378,137.0 m, which corresponds + to the WGS84 ellipsoid) + flattening : float + The flattening of the ellipsoid used to represent the Earth, usually + between 1/250 and 1/150 (default is 1/298.257223563, which + corresponds to the WGS84 ellipsoid) Returns ------- radius : float Earth Radius at the desired latitude in meters """ - # Select the desired datum (i.e. the ellipsoid parameters) - flattening = self.ellipsoid.flattening - semi_major_axis = self.ellipsoid.semi_major_axis - - # Calculate the semi minor axis length - # semi_minor_axis = semi_major_axis - semi_major_axis*(flattening**(-1)) semi_minor_axis = semi_major_axis * (1 - flattening) - # Convert latitude to radians + # Numpy trigonometric functions work with radians, so convert to radians lat = lat * np.pi / 180 # Calculate the Earth Radius in meters @@ -3555,9 +3584,6 @@ def calculate_earth_radius(self, lat): ) ) - # Convert latitude to degrees - lat = lat * 180 / np.pi - return e_radius @staticmethod diff --git a/tests/test_environment.py b/tests/test_environment.py index 818e43303..ddf19ef7d 100644 --- a/tests/test_environment.py +++ b/tests/test_environment.py @@ -371,20 +371,57 @@ def test_export_environment(example_env_robust): os.remove("environment.json") -def test_utm_to_geodesic(example_env_robust): - """Tests the conversion from UTM to geodesic coordinates. +def test_geodesic_to_utm(): + """Tests the conversion from geodesic to UTM coordinates.""" + x, y, utm_zone, utm_letter, hemis, EW = Environment.geodesic_to_utm( + lat=32.990254, + lon=-106.974998, + semi_major_axis=6378137.0, # WGS84 + flattening=1 / 298.257223563, # WGS84 + ) + assert np.isclose(x, 315468.64, atol=1e-5) == True + assert np.isclose(y, 3651938.65, atol=1e-5) == True + assert utm_zone == 13 + assert utm_letter == "S" + assert hemis == "N" + assert EW == "W" + + +def test_utm_to_geodesic(): + """Tests the conversion from UTM to geodesic coordinates.""" + + lat, lon = Environment.utm_to_geodesic( + x=315468.64, + y=3651938.65, + utm_zone=13, + hemis="N", + semi_major_axis=6378137.0, # WGS84 + flattening=1 / 298.257223563, # WGS84 + ) + assert np.isclose(lat, 32.99025, atol=1e-5) == True + assert np.isclose(lon, -106.9750, atol=1e-5) == True + + +@pytest.mark.parametrize( + "lat, radius", [(0, 6378137.0), (90, 6356752.31424518), (-90, 6356752.31424518)] +) +def test_earth_radius_calculation(lat, radius): + """Tests if the earth radius calculation is correct. It takes 3 different + latitudes and their expected results and compares them with the results + from the method. Parameters ---------- - example_env_robust : rocketpy.Environment - Example environment object to be tested. + lat : float + The latitude in decimal degrees. + radius : float + The expected radius in meters at the given latitude. """ - # TODO: maybe the utm_to_geodesic should be a static method, or move to tools.py - lat, lon = example_env_robust.utm_to_geodesic( - x=315468.64, y=3651938.65, utm_zone=13, hemis="N" - ) - assert np.isclose(lat, 32.99025, atol=1e-5) == True - assert np.isclose(lon, -106.9750, atol=1e-5) == True + semi_major_axis = 6378137.0 # WGS84 + flattening = 1 / 298.257223563 # WGS84 + res = Environment.calculate_earth_radius(lat, semi_major_axis, flattening) + assert pytest.approx(res, abs=1e-8) == radius + @pytest.mark.parametrize( "angle, deg, arc_min, arc_sec", From 316ee46c8d3e40c5b94b4c5c90bd551150fe5a23 Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Fri, 22 Sep 2023 17:39:34 -0300 Subject: [PATCH 4/4] DEP: remove add_fin method --- rocketpy/Rocket.py | 39 --------------------------------------- 1 file changed, 39 deletions(-) diff --git a/rocketpy/Rocket.py b/rocketpy/Rocket.py index 35ed231f1..9f6b1b44c 100644 --- a/rocketpy/Rocket.py +++ b/rocketpy/Rocket.py @@ -1199,42 +1199,3 @@ def all_info(self): self.plots.all() return None - - def add_fin( - self, - number_of_fins=4, - cl=2 * np.pi, - cpr=1, - cpz=1, - gammas=[0, 0, 0, 0], - angular_positions=None, - ): - self.aerodynamic_surfaces = Components() - pi = np.pi - # Calculate angular positions if not given - if angular_positions is None: - angular_positions = ( - np.array(range(number_of_fins)) * 2 * pi / number_of_fins - ) - else: - angular_positions = np.array(angular_positions) * pi / 180 - # Convert gammas to degree - if isinstance(gammas, (int, float)): - gammas = [(pi / 180) * gammas for i in range(number_of_fins)] - else: - gammas = [(pi / 180) * gamma for gamma in gammas] - for i in range(number_of_fins): - # Get angular position and inclination for current fin - angular_position = angular_positions[i] - gamma = gammas[i] - # Calculate position vector - cpx = cpr * np.cos(angular_position) - cpy = cpr * np.sin(angular_position) - position_vector = np.array([cpx, cpy, cpz]) - # Calculate chord vector - aux_vector = np.array([cpy, -cpx, 0]) / (cpr) - chord_vector = ( - np.cos(gamma) * np.array([0, 0, 1]) - np.sin(gamma) * aux_vector - ) - self.aerodynamic_surfaces.append([position_vector, chord_vector]) - return None