diff --git a/src/orbit_visualiser/core/__init__.py b/src/orbit_visualiser/core/__init__.py index eb1bcd0..9b0a054 100644 --- a/src/orbit_visualiser/core/__init__.py +++ b/src/orbit_visualiser/core/__init__.py @@ -4,11 +4,11 @@ from .astrodynamics.types import OrbitType from .astrodynamics.keplerian.classification import orbit_type -from .astrodynamics.keplerian.anomlies import mean_anomaly, eccentric_anomaly +from .astrodynamics.keplerian.anomalies import mean_anomaly, eccentric_anomaly from .astrodynamics.keplerian.elements import (eccentricity_vector_from_state, eccentricity_from_state, true_anomaly_from_state, semi_parameter_from_momentum, semi_parameter_from_eccentricity, semimajor_axis, - semiminor_axis, periapsis, apoapsis, + semiminor_axis, radius_of_periapsis, radius_of_apoapsis, asymptote_anomaly, turning_angle, aiming_radius, orbital_period, mean_motion) from .astrodynamics.keplerian.dynamics import (specific_ang_momentum_from_state, specific_ang_momentum, diff --git a/src/orbit_visualiser/core/astrodynamics/keplerian/anomlies.py b/src/orbit_visualiser/core/astrodynamics/keplerian/anomalies.py similarity index 90% rename from src/orbit_visualiser/core/astrodynamics/keplerian/anomlies.py rename to src/orbit_visualiser/core/astrodynamics/keplerian/anomalies.py index 9075ae9..655b6a3 100644 --- a/src/orbit_visualiser/core/astrodynamics/keplerian/anomlies.py +++ b/src/orbit_visualiser/core/astrodynamics/keplerian/anomalies.py @@ -12,12 +12,12 @@ def eccentric_anomaly(e: float, nu: float) -> float: e : float Eccentricity nu : float - True anomaly (rads) + True anomaly (rad) Returns ------- float - The eccentric anomaly (rads) + The eccentric anomaly (rad) """ type = orbit_type(e) if type is OrbitType.CIRCULAR: @@ -45,15 +45,15 @@ def mean_anomaly(e: float, nu: float, e_anomaly: float) -> float: e : float Eccentricity nu : float - True anomaly (rads) + True anomaly (rad) e_anomaly : float - The eccentric anomaly (rads) + The eccentric anomaly (rad) Returns ------- float - The mean anomaly (rads) + The mean anomaly (rad) """ type = orbit_type(e) if type in (OrbitType.CIRCULAR, OrbitType.ELLIPTICAL): diff --git a/src/orbit_visualiser/core/astrodynamics/keplerian/elements.py b/src/orbit_visualiser/core/astrodynamics/keplerian/elements.py index b88695d..01feede 100644 --- a/src/orbit_visualiser/core/astrodynamics/keplerian/elements.py +++ b/src/orbit_visualiser/core/astrodynamics/keplerian/elements.py @@ -148,7 +148,7 @@ def semiminor_axis(e: float, a: float) -> float: return np.nan -def periapsis(p: float, e: float) -> float: +def radius_of_periapsis(p: float, e: float) -> float: """ Calculates the radius of periapsis from the semi-parameter and eccentricity. @@ -166,7 +166,26 @@ def periapsis(p: float, e: float) -> float: """ return p/(1 + e) -def apoapsis(e: float, a: float) -> float: +def periapsis(p: float, e: NDArray[np.float64]) -> NDArray[np.float64]: + """ + Calculate the periapsis vector from the eccentricity vector and the semi-parameter. + + Parameters + ---------- + p : float + The semi-parameter (km) + e : NDArray[np.float64] + The eccentricity vector + + Returns + ------- + NDArray[np.float64] + Periapsis vector (km) + """ + e_norm = np.linalg.norm(e) + return p/((e_norm(1 + e_norm)))*e + +def radius_of_apoapsis(e: float, a: float) -> float: """ Calculates the radius of apoapsis using the eccentricity and semi-major axis. @@ -187,6 +206,29 @@ def apoapsis(e: float, a: float) -> float: return a*(1 + e) +def apoapsis(e: NDArray[np.float64], a: float) -> NDArray[np.float64]: + """ + Calculate the apoapsis vector from the eccentricity vector and the semi-major axis. + + Parameters + ---------- + e : NDArray[np.float64] + Eccentricity vector + a : float + Semi-major axis (km) + + Returns + ------- + NDArray[np.float64] + Apoapsis vector (km) + """ + + e_norm = np.linalg.norm(e) + if orbit_type(e_norm) is OrbitType.PARABOLIC: + return np.empty(e.shape) + + return -(a*(1 + e_norm)/e_norm)*e + def asymptote_anomaly(e: float) -> float: """ Calculate the asymptote of the true anomaly for open orbits using the eccentricity. diff --git a/src/orbit_visualiser/core/astrodynamics/keplerian/state.py b/src/orbit_visualiser/core/astrodynamics/keplerian/state.py index 488536d..ebe6533 100644 --- a/src/orbit_visualiser/core/astrodynamics/keplerian/state.py +++ b/src/orbit_visualiser/core/astrodynamics/keplerian/state.py @@ -1,7 +1,7 @@ from math import pi import numpy as np from numpy.typing import NDArray -from typing import Callable, Literal +from typing import Literal from orbit_visualiser.core.astrodynamics.types import OrbitType from orbit_visualiser.core.astrodynamics.keplerian.elements import semi_parameter_from_eccentricity from orbit_visualiser.core.astrodynamics.keplerian.dynamics import specific_ang_momentum @@ -49,7 +49,7 @@ def state_pf_from_e_rp( elif state == "both": return [r, v] -def perifocal_position(e: float, p: float, nu: float) -> NDArray[np.float64]: +def perifocal_position(e: float, p: float, nu: float | NDArray[np.float64]) -> NDArray[np.float64]: """ Perifocal orbit equation. @@ -59,17 +59,18 @@ def perifocal_position(e: float, p: float, nu: float) -> NDArray[np.float64]: Eccentricity p : float Semi-parameter (km) - nu : float - True anomaly (rads) + nu : float | NDArray[np.float64] + True anomaly (rad) Returns ------- NDArray[np.float64] The numpy array of the perifocal orbital position """ - return p*(1/(1 + e*np.cos(nu)))*np.array([np.cos(nu), np.sin(nu)]) + z = 0 if isinstance(nu, float) else np.zeros(nu.shape) + return p*(1/(1 + e*np.cos(nu)))*np.array([np.cos(nu), np.sin(nu), z]) -def perifocal_velocity(e: float, mu: float, h: float, nu: float) -> NDArray[np.float64]: +def perifocal_velocity(e: float, mu: float, h: float, nu: float | NDArray[np.float64]) -> NDArray[np.float64]: """ The perifocal velocity equation. @@ -81,15 +82,16 @@ def perifocal_velocity(e: float, mu: float, h: float, nu: float) -> NDArray[np.f Gravitational parameter (km^3/s^2) h : float Specific angular momentum (km^2/s) - nu : float - True anomaly (rads) + nu : float | NDArray[np.float64] + True anomaly (rad) Returns ------- NDArray[np.float64] The numpy array of the perifocal velocity """ - return (mu/h)*np.array([-np.sin(nu), e + np.cos(nu)]) + z = 0 if isinstance(nu, float) else np.zeros(nu.shape) + return (mu/h)*np.array([-np.sin(nu), e + np.cos(nu), z]) def radial_azimuthal_velocity( @@ -106,7 +108,7 @@ def radial_azimuthal_velocity( Parameters ---------- nu : float - True anomaly (rads) + True anomaly (rad) mu : float Gravitational parameter (km^3/s^2) h : float @@ -114,7 +116,7 @@ def radial_azimuthal_velocity( e : float Eccentricity asymp_anomaly : float - The asymptote of the free anomaly (rads) + The asymptote of the free anomaly (rad) Returns ------- @@ -166,9 +168,9 @@ def radius_from_orbit_eq(nu: float, asymp_anomaly: float, p: float, e: float) -> Parameters ---------- nu : float - True anomaly (rads) + True anomaly (rad) asymp_anomaly : float - The asymptote of the free anomaly (rads) + The asymptote of the free anomaly (rad) p : float Semi-parameter (km) e : float @@ -191,9 +193,9 @@ def escape_velocity(nu: float, asymp_anomaly: float, mu: float, r: float) -> flo Parameters ---------- nu : float - The true anomaly (rads) + The true anomaly (rad) asymp_anomaly : float - The true anomaly of the asymptot (rads) + The true anomaly of the asymptote (rad) mu : float Gravitational parameter (km^3/s^2) r : float @@ -216,16 +218,16 @@ def flight_angle(nu: float, asymp_anomaly: float, e: float) -> float: Parameters ---------- nu : float - True anomaly (rads) + True anomaly (rad) asymp_anomaly : float - The true anomaly of the asymptot (rads) + The true anomaly of the asymptot (rad) e : float Eccentricity Returns ------- float - The satellite flight angle (rads) + The satellite flight angle (rad) """ if np.isclose(abs(nu), asymp_anomaly): return np.sign(nu)*pi/2 @@ -239,7 +241,7 @@ def time_since_periapsis(m_anomaly: float, period: float, p: float, h: float, e: Parameters ---------- m_anomaly : float - The mean anomaly (rads) + The mean anomaly (rad) period : float The orbital period (s) p : float diff --git a/src/orbit_visualiser/core/orbit.py b/src/orbit_visualiser/core/orbit.py index 034ec5d..d1ccd43 100644 --- a/src/orbit_visualiser/core/orbit.py +++ b/src/orbit_visualiser/core/orbit.py @@ -5,8 +5,8 @@ from numpy.typing import NDArray from orbit_visualiser.core.astrodynamics.keplerian.state import state_pf_from_e_rp from orbit_visualiser.core.astrodynamics.keplerian.elements import (eccentricity_from_state, semi_parameter_from_momentum, - periapsis, semimajor_axis, semiminor_axis, - apoapsis, asymptote_anomaly, turning_angle, + radius_of_periapsis, semimajor_axis, semiminor_axis, + radius_of_apoapsis, asymptote_anomaly, turning_angle, aiming_radius, orbital_period, mean_motion) from orbit_visualiser.core.astrodynamics.keplerian.dynamics import (specific_orbital_energy, characteristic_energy, excess_velocity, specific_ang_momentum_from_state) @@ -58,7 +58,7 @@ class Orbit(): mu : float The gravitational parameter of the central body (km^3/s^2) nu : float - The true anomaly of the satellite (rads) + The true anomaly of the satellite (rad) """ position : Sequence | NDArray[np.float64] velocity : Sequence | NDArray[np.float64] @@ -78,11 +78,11 @@ def orbit_type(self) -> OrbitType: @cached_property def semi_parameter(self) -> float: - return semi_parameter_from_momentum(specific_ang_momentum_from_state(self.position, self.velocity), self.mu) + return semi_parameter_from_momentum(np.linalg.norm(specific_ang_momentum_from_state(self.position, self.velocity)), self.mu) @cached_property def radius_of_periapsis(self) -> float: - return periapsis(self.semi_parameter, self.eccentricity) + return radius_of_periapsis(self.semi_parameter, self.eccentricity) @cached_property def semimajor_axis(self) -> float: @@ -94,7 +94,7 @@ def semiminor_axis(self) -> float: @cached_property def radius_of_apoapsis(self) -> float: - return apoapsis(self.eccentricity, self.semimajor_axis) + return radius_of_apoapsis(self.eccentricity, self.semimajor_axis) @cached_property def asymptote_anomaly(self) -> float: diff --git a/src/orbit_visualiser/core/propagation.py b/src/orbit_visualiser/core/propagation.py index ce7775d..4eda378 100644 --- a/src/orbit_visualiser/core/propagation.py +++ b/src/orbit_visualiser/core/propagation.py @@ -27,25 +27,23 @@ def two_body_pf_ode(mu: float, t: float, state: NDArray[np.float64], ) -> NDArra NDArray[np.float64] The values to be numerically integrated [dx, d^2x] """ - x, y, v_x, v_y = state + pos = state[:3] + r = np.linalg.norm(pos) - r = np.hypot(x, y) + vel = state[3:] + acc: NDArray[np.float64] = -(mu/r**3)*pos - a_x = -(mu/r**3)*x - a_y = -(mu/r**3)*y - - return np.array([v_x, v_y, a_x, a_y]) + return np.concatenate((vel, acc)) def get_init_conditions_from_orbit(orbit: Orbit) -> NDArray[np.float64]: """ - Takes orbital elements and returns the initial conditions (position and velocity) for orbit propagation. + Takes orbit object and returns the initial conditions (position and velocity) for orbit propagation. Parameters ---------- orbit: NewOrbit Orbit object to get the initial conditions from - Returns ------- NDArray[np.float64] diff --git a/src/orbit_visualiser/core/satellite.py b/src/orbit_visualiser/core/satellite.py index 36b0c26..473a811 100644 --- a/src/orbit_visualiser/core/satellite.py +++ b/src/orbit_visualiser/core/satellite.py @@ -3,7 +3,7 @@ from typing import Sequence from orbit_visualiser.core.astrodynamics.keplerian.state import (speed, radial_azimuthal_velocity, escape_velocity, radius_from_state, flight_angle, time_since_periapsis) -from orbit_visualiser.core.astrodynamics.keplerian.anomlies import mean_anomaly, eccentric_anomaly +from orbit_visualiser.core.astrodynamics.keplerian.anomalies import mean_anomaly, eccentric_anomaly from orbit_visualiser.core.astrodynamics.keplerian.dynamics import specific_ang_momentum_from_state from orbit_visualiser.core.astrodynamics.keplerian.elements import true_anomaly_from_state from orbit_visualiser.core.orbit import CentralBody, Orbit @@ -93,7 +93,7 @@ def radial_azimuthal_velocity(self) -> NDArray[np.float64]: return radial_azimuthal_velocity( self.true_anomaly, self.central_body.mu, - self.specific_angular_momentum, + np.linalg.norm(self.specific_angular_momentum), orbit.eccentricity, orbit.asymptote_anomaly ) @@ -141,6 +141,6 @@ def time_since_periapsis(self) -> float: self.mean_anomaly, orbit.orbital_period, orbit.semi_parameter, - self.specific_angular_momentum, + np.linalg.norm(self.specific_angular_momentum), orbit.eccentricity ) \ No newline at end of file diff --git a/src/orbit_visualiser/ui/config/config_controller.py b/src/orbit_visualiser/ui/config/config_controller.py index 704bd32..3b31395 100644 --- a/src/orbit_visualiser/ui/config/config_controller.py +++ b/src/orbit_visualiser/ui/config/config_controller.py @@ -1,7 +1,7 @@ from tkinter import Event from typing import Callable from orbit_visualiser.core import Orbit, Satellite, CentralBody -from orbit_visualiser.ui.config.config_builder import OrbitConfigBuilder, OrbitConfigBuilder +from orbit_visualiser.ui.config.config_builder import OrbitConfigBuilder from orbit_visualiser.ui.config.variables_panel.variables_panel_controller import VariablesController from orbit_visualiser.ui.config.properties_panel.properties_panel_controller import PropertiesController #from orbit_visualiser.ui.config.display_panel.display_panel_controller import DisplayController diff --git a/src/orbit_visualiser/ui/config/properties_panel/properties_panel_builder.py b/src/orbit_visualiser/ui/config/properties_panel/properties_panel_builder.py index 7d7bdc5..176a119 100644 --- a/src/orbit_visualiser/ui/config/properties_panel/properties_panel_builder.py +++ b/src/orbit_visualiser/ui/config/properties_panel/properties_panel_builder.py @@ -40,7 +40,7 @@ def __init__( "e_anomaly" : PropertySpec("Eccentric anomaly", "°", lambda sat: np.degrees(sat.eccentric_anomaly)), "m_anomaly" : PropertySpec("Mean anomaly", "°", lambda sat: np.degrees(sat.mean_anomaly)), "time_periapsis" : PropertySpec("Time since periapsis", "s", lambda sat: sat.time_since_periapsis), - "ang_momentum" : PropertySpec("Angular momentum", "km²/s", lambda sat: sat.specific_angular_momentum), + "ang_momentum" : PropertySpec("Angular momentum", "km²/s", lambda sat: np.linalg.norm(sat.specific_angular_momentum)), "speed" : PropertySpec("Orbital speed", "km/s", lambda sat: sat.speed), "azim_velocity" : PropertySpec("Azimuthal velocity", "km/s", lambda sat: sat.radial_azimuthal_velocity[1]), "radial_velocity" : PropertySpec("Radial velocity", "km/s", lambda sat: sat.radial_azimuthal_velocity[0]), diff --git a/src/orbit_visualiser/ui/figure/orbit_figure.py b/src/orbit_visualiser/ui/figure/orbit_figure.py index a9baa6f..2735a84 100644 --- a/src/orbit_visualiser/ui/figure/orbit_figure.py +++ b/src/orbit_visualiser/ui/figure/orbit_figure.py @@ -7,7 +7,6 @@ import numpy as np from numpy.typing import NDArray from orbit_visualiser.core import Orbit, Satellite, OrbitType, perifocal_position -from orbit_visualiser.core.astrodynamics.types import OrbitType # TODO: Fix bug where scroll zoom doesn't register as changing the view so the native matplotlib home button has unexpected (and often undesirable) behaviour. class OrbitFigure(): @@ -78,7 +77,7 @@ def _initialise_plot(self) -> None: # Plot the initial orbit orbit = self._satellite.orbit nu = self._get_anomaly_data(orbit) - x, y = perifocal_position(orbit.eccentricity, orbit.semi_parameter, nu) + x, y, _ = perifocal_position(orbit.eccentricity, orbit.semi_parameter, nu) self._line, = self._ax.plot(x, y, color = "#2F2F2F", alpha = 0.5, linewidth = 1.5) # Plot the central body @@ -94,7 +93,7 @@ def _initialise_plot(self) -> None: #self.plot_periapsis_point() - f = self._zoom_factory(self._ax, 1.1) + self._zoom_factory(self._ax, 1.1) def _build_canvas(self) -> None: self._canvas = FigureCanvasTkAgg(self._fig, master = self._figure_frame) @@ -123,13 +122,13 @@ def _get_anomaly_data(self, orbit: Orbit) -> NDArray[np.float64]: def redraw_orbit(self) -> None: orbit = self._satellite.orbit nu = self._get_anomaly_data(orbit) - x, y = perifocal_position(orbit.eccentricity, orbit.semi_parameter, nu) + x, y, _ = perifocal_position(orbit.eccentricity, orbit.semi_parameter, nu) self._line.set_data(x, y) self._canvas.draw_idle() def redraw_satellite(self) -> None: - x, y = self._satellite.position + x, y, _ = self._satellite.position self._sat_point.set_data((x,), (y,)) self._canvas.draw_idle() diff --git a/tests/propagation/propagation_test.py b/tests/propagation/propagation_test.py index 581a62e..0fcbfbb 100644 --- a/tests/propagation/propagation_test.py +++ b/tests/propagation/propagation_test.py @@ -47,7 +47,7 @@ def test_1p_propagation_pos( pos_at_rp = init_conditions[:2] prop_pos_at_rp = prop.y[:2, -1] - assert np.allclose(pos_at_rp, prop_pos_at_rp, atol = 0.01) + assert np.allclose(pos_at_rp, prop_pos_at_rp, atol = 0.05) @pytest.mark.parametrize("e, rp, mu", typical_closed_test_cases) @@ -62,11 +62,10 @@ def test_1p_phase_shift( typical closed orbits is within an acceptable tolerance to the analytical solution. """ orbit: Orbit = orbit_factory_from_elements(e, rp, mu, 0.0) - init_conditions = get_init_conditions_from_orbit(orbit) prop = run_orbit_prop(orbit, orbit.orbital_period) prop_pos_at_rp = prop.y[:2, -1] phase_shift = np.arctan2(prop_pos_at_rp[1], prop_pos_at_rp[0]) - assert np.isclose(phase_shift, 0, atol=1e-8) \ No newline at end of file + assert np.isclose(phase_shift, 0, atol=5e-8) \ No newline at end of file diff --git a/tests/satellite_test.py b/tests/satellite_test.py index 8af7bb5..2d6bed5 100644 --- a/tests/satellite_test.py +++ b/tests/satellite_test.py @@ -53,7 +53,7 @@ def test_satellite_azimuthal_velocity_sanity( satellite: Satellite = satellite_factory_from_elements(e = e, rp = rp, mu = mu, nu = nu) v_azim_vals = [ - satellite.specific_angular_momentum/satellite.radius, + np.linalg.norm(satellite.specific_angular_momentum)/satellite.radius, satellite.speed*np.cos(satellite.flight_angle) ] diff --git a/tests/unit_tests/dynamics_test.py b/tests/unit_tests/dynamics_test.py index 44d013d..2920952 100644 --- a/tests/unit_tests/dynamics_test.py +++ b/tests/unit_tests/dynamics_test.py @@ -6,9 +6,9 @@ vis_viva_speed) @pytest.mark.parametrize("r, v, expected", [ - (np.array([1.0, 0]), np.array([0, 1.0]), np.array(1.0)), - (np.array([1.0, 2.0]), np.array([2.0, 4.0]), np.array(0.0)), - (np.array([50_000, 0]), np.array([0, 4.0]), np.array(4*50_000)) + (np.array([1.0, 0, 0]), np.array([0, 1.0, 0]), np.array([0.0, 0.0, 1.0])), + (np.array([1.0, 2.0, 0]), np.array([2.0, 4.0, 0]), np.array([0.0, 0.0, 0.0])), + (np.array([50_000, 0, 0]), np.array([0, 4.0, 0]), np.array([0.0, 0.0, 4*50_000.0])) ]) def test_specific_angular_momentum_from_state( r: NDArray[np.float64], diff --git a/tests/unit_tests/elements_test.py b/tests/unit_tests/elements_test.py index f6244f9..98282da 100644 --- a/tests/unit_tests/elements_test.py +++ b/tests/unit_tests/elements_test.py @@ -5,7 +5,7 @@ from orbit_visualiser.core import (OrbitType, eccentricity_vector_from_state, eccentricity_from_state, true_anomaly_from_state, semi_parameter_from_momentum, semi_parameter_from_eccentricity, semimajor_axis, - semiminor_axis, periapsis, apoapsis, + semiminor_axis, radius_of_periapsis, radius_of_apoapsis, asymptote_anomaly, turning_angle, aiming_radius, orbital_period, mean_motion) @@ -111,7 +111,7 @@ def test_periapsis(p: float, e: float, expected: float): Test that the formula for the periapsis using the eccentricity and semi-parameter gives the expected value. """ - result = periapsis(p, e) + result = radius_of_periapsis(p, e) assert np.isclose(result, expected) @pytest.mark.parametrize("e, a, expected", [ @@ -123,7 +123,7 @@ def test_apoapsis(e: float, a: float, expected: float): Test that the formula for the periapsis using the eccentricity and semi-major axis gives the expected value. """ - result = apoapsis(e, a) + result = radius_of_apoapsis(e, a) assert np.isclose(result, expected, equal_nan = True) @pytest.mark.parametrize("e, expected", [ diff --git a/tests/unit_tests/state_test.py b/tests/unit_tests/state_test.py index b5ada90..2600fe4 100644 --- a/tests/unit_tests/state_test.py +++ b/tests/unit_tests/state_test.py @@ -7,7 +7,7 @@ flight_angle, time_since_periapsis) @pytest.mark.parametrize("e, p, nu, expected", [ - (0.5, 75_000.0, pi/2, np.array([0.0, 75_000.0])) + (0.5, 75_000.0, pi/2, np.array([0.0, 75_000.0, 0])) ]) def test_perifocal_position(e: float, p: float, nu: float, expected: NDArray[np.float64]): """ @@ -17,7 +17,7 @@ def test_perifocal_position(e: float, p: float, nu: float, expected: NDArray[np. assert np.allclose(result, expected) @pytest.mark.parametrize("e, mu, h, nu, expected", [ - (0.5, 398_600.0, 199_300, pi/6, np.array([-1.0, 1 + np.sqrt(3)])) + (0.5, 398_600.0, 199_300, pi/6, np.array([-1.0, 1 + np.sqrt(3), 0])) ]) def test_perifocal_velocity(e: float, mu: float, h:float, nu: float, expected: NDArray[np.float64]): """