diff --git a/crates/lox-orbits/src/python.rs b/crates/lox-orbits/src/python.rs index 992a1943..7e4c5a3c 100644 --- a/crates/lox-orbits/src/python.rs +++ b/crates/lox-orbits/src/python.rs @@ -41,6 +41,7 @@ use crate::propagators::semi_analytical::{Vallado, ValladoError}; use crate::propagators::sgp4::{Sgp4, Sgp4Error}; use crate::propagators::Propagator; use crate::states::ToCartesian; +use crate::trajectories::TrajectoryTransformationError; use crate::{ frames::FrameTransformationProvider, states::State, @@ -49,6 +50,13 @@ use crate::{ mod generated; +impl From for PyErr { + fn from(err: TrajectoryTransformationError) -> Self { + // FIXME: wrong error type + PyValueError::new_err(err.to_string()) + } +} + impl From for PyErr { fn from(err: FindEventError) -> Self { // FIXME: wrong error type @@ -382,9 +390,19 @@ impl PyState { } fn to_origin(&self, target: &Bound<'_, PyAny>, ephemeris: &Bound<'_, PySpk>) -> PyResult { + if self.0.reference_frame() != PyFrame::Icrf { + return Err(PyValueError::new_err( + "only inertial frames are supported for conversion to Keplerian elements", + )); + } let target = PyBody::try_from(target)?; - let spk = ephemeris.borrow(); - todo!() + let spk = &ephemeris.borrow().0; + let s1 = self + .0 + .with_frame(Icrf) + .to_origin(target, spk)? + .with_frame(PyFrame::Icrf); + Ok(Self(s1)) } fn to_keplerian(&self) -> PyResult { @@ -620,6 +638,32 @@ impl PyTrajectory { } Err(PyValueError::new_err("invalid time argument")) } + + #[pyo3(signature = (frame, provider=None))] + fn to_frame( + &self, + frame: PyFrame, + provider: Option<&Bound<'_, PyUt1Provider>>, + ) -> PyResult { + let mut states: Vec> = + Vec::with_capacity(self.0.states().len()); + for s in self.0.states() { + states.push(PyState(s).to_frame(frame.clone(), provider)?.0); + } + Ok(PyTrajectory(Trajectory::new(&states)?)) + } + + fn to_origin(&self, target: &Bound<'_, PyAny>, ephemeris: &Bound<'_, PySpk>) -> PyResult { + let target = PyBody::try_from(target)?; + let spk = &ephemeris.borrow().0; + let s1 = self + .0 + .clone() + .with_frame(Icrf) + .to_origin(target, spk)? + .with_frame(PyFrame::Icrf); + Ok(Self(s1)) + } } #[pyclass(name = "Event", module = "lox_space", frozen)] diff --git a/crates/lox-space/lox_space.pyi b/crates/lox-space/lox_space.pyi index 5d3257bd..bbc5242e 100644 --- a/crates/lox-space/lox_space.pyi +++ b/crates/lox-space/lox_space.pyi @@ -7,498 +7,324 @@ type Epoch = Literal["jd", "mjd", "j1950", "j2000"] type Unit = Literal["seconds", "days", "centuries"] type Vec3 = tuple[float, float, float] - def find_events(func: Callable[[float], float], start: Time, times: list[float]): ... - - def find_windows( - func: Callable[[float], float], start: Time, end: Time, times: list[float] + func: Callable[[float], float], start: Time, end: Time, times: list[float] ): ... - - def visibility( - times: list[Time], - gs: GroundLocation, - min_elevation: float, - sc: Trajectory, - provider: UT1Provider, + times: list[Time], + gs: GroundLocation, + min_elevation: float, + sc: Trajectory, + provider: UT1Provider, ): ... - class Frame: def __new__(cls, abbreviation: str): ... - def name(self) -> str: ... - def abbreviation(self) -> str: ... +class SPK: + def __new__(cls, path): ... class State: def __new__( - cls, - time: Time, - position: Vec3, - velocity: Vec3, - origin: Origin | None = None, - frame: Frame | None = None, + cls, + time: Time, + position: Vec3, + velocity: Vec3, + origin: Origin | None = None, + frame: Frame | None = None, ): ... - def time(self) -> Time: ... - def origin(self) -> Origin: ... - def reference_frame(self) -> Frame: ... - def position(self) -> np.ndarray: ... - def velocity(self) -> np.ndarray: ... - - def to_frame(self, frame: Frame) -> State: ... - + def to_frame(self, frame: Frame) -> Self: ... + def to_origin(self, target: Origin, ephemeris: SPK) -> Self: ... def to_keplerian(self) -> Keplerian: ... - def rotation_lvlh(self) -> np.ndarray: ... - def to_ground_location(self) -> GroundLocation: ... - class Keplerian: def __new__( - cls, - time: Time, - semi_major_axis: float, - eccentricity: float, - inclination: float, - longitude_of_ascending_node: float, - argument_of_periapsis: float, - true_anomaly: float, - origin: Origin | None = None, + cls, + time: Time, + semi_major_axis: float, + eccentricity: float, + inclination: float, + longitude_of_ascending_node: float, + argument_of_periapsis: float, + true_anomaly: float, + origin: Origin | None = None, ): ... - def time(self) -> Time: ... - def origin(self) -> Origin: ... - def semi_major_axis(self) -> float: ... - def eccentricity(self) -> float: ... - def inclination(self) -> float: ... - def longitude_of_ascending_node(self) -> float: ... - def argument_of_periapsis(self) -> float: ... - def true_anomaly(self) -> float: ... - def to_cartesian(self) -> State: ... - def orbital_period(self) -> TimeDelta: ... - class Trajectory: def __new__(cls, states: list[State]): ... - @classmethod def from_numpy( - cls, - start_time: Time, - states: np.ndarray, - origin: Origin | None = None, - frame: Frame | None = None, + cls, + start_time: Time, + states: np.ndarray, + origin: Origin | None = None, + frame: Frame | None = None, ) -> Self: ... - def origin(self) -> Origin: ... - def reference_frame(self) -> Frame: ... - def to_numpy(self) -> np.ndarray: ... - def states(self) -> list[State]: ... - def find_events(self, func: Callable[[State], float]) -> list[Event]: ... - def find_windows(self, func: Callable[[State], float]) -> list[Window]: ... - def interpolate(self, time: Time) -> State: ... - + def to_frame(self, frame: Frame) -> Self: ... + def to_origin(self, target: Origin, ephemeris: SPK) -> Self: ... class Event: def time(self) -> Time: ... - def cross(self) -> str: ... - class Window: def start(self) -> Time: ... - def end(self) -> Time: ... - def duration(self) -> TimeDelta: ... - class Vallado: def __new__(cls, initial_state: State, max_iter: int | None = None): ... - @overload def propagate(self, time: Time) -> State: ... - @overload def propagate(self, time: list[Time]) -> Trajectory: ... - def propagate(self, time: Time | list[Time]) -> State | Trajectory: ... - class GroundLocation: def __new__( - cls, - planet: Planet, - longitude: float, - latitude: float, - altitude: float, + cls, + planet: Planet, + longitude: float, + latitude: float, + altitude: float, ): ... - def longitude(self) -> float: ... - def latitude(self) -> float: ... - def altitude(self) -> float: ... - def observables(self) -> Observables: ... - def rotation_to_topocentric(self) -> np.ndarray: ... - class GroundPropagator: def __new__(cls, location: GroundLocation, provider: UT1Provider): ... - @overload def propagate(self, time: Time) -> State: ... - @overload def propagate(self, time: list[Time]) -> Trajectory: ... - def propagate(self, time: Time | list[Time]) -> State | Trajectory: ... - class SGP4: def __new__(cls, tle: str): ... - def time(self) -> Time: ... - @overload def propagate(self, time: Time) -> State: ... - @overload def propagate(self, time: list[Time]) -> Trajectory: ... - def propagate(self, time: Time | list[Time]) -> State | Trajectory: ... - class Observables: - def __new__(cls, azimuth: float, elevation: float, range: float, range_rate: float): ... - + def __new__( + cls, azimuth: float, elevation: float, range: float, range_rate: float + ): ... def azimuth(self) -> float: ... - def elevation(self) -> float: ... - def range(self) -> float: ... - def range_rate(self) -> float: ... - class Time: def __new__( - cls, - scale: Scale, - year: int, - month: int, - day: int, - hour: int = 0, - minute: int = 0, - seconds: float = 0.0, + cls, + scale: Scale, + year: int, + month: int, + day: int, + hour: int = 0, + minute: int = 0, + seconds: float = 0.0, ): ... - @classmethod def from_julian_date(cls, scale: Scale, jd: float, epoch: str = "jd") -> Self: ... - @classmethod def from_two_part_julian_date( - cls, scale: Scale, jd1: float, jd2: float + cls, scale: Scale, jd1: float, jd2: float ) -> Self: ... - @classmethod def from_day_of_year( - cls, - scale: Scale, - year: int, - doy: int, - hour: int = 0, - minute: int = 0, - seconds: float = 0.0, + cls, + scale: Scale, + year: int, + doy: int, + hour: int = 0, + minute: int = 0, + seconds: float = 0.0, ) -> Self: ... - @classmethod def from_iso(cls, iso: str, scale: Scale | None = None) -> Self: ... - @classmethod def from_seconds(cls, scale: Scale, seconds: int, subsecond: float) -> Self: ... - def seconds(self) -> int: ... - def subsecond(self) -> float: ... - def __str__(self) -> str: ... - def __repr__(self) -> str: ... - def __add__(self, other: TimeDelta) -> Self: ... - @overload def __sub__(self, other: TimeDelta) -> Self: ... - @overload def __sub__(self, other: Time) -> TimeDelta: ... - def __eq__(self, other: object) -> bool: ... - def __lt__(self, other: object) -> bool: ... - def __le__(self, other: object) -> bool: ... - def isclose( - self, other: Time, rel_tol: float = 1e-8, abs_tol: float = 1e-14 + self, other: Time, rel_tol: float = 1e-8, abs_tol: float = 1e-14 ) -> bool: ... - def julian_date( - self, - epoch: Epoch = "jd", - unit: Unit = "days", + self, + epoch: Epoch = "jd", + unit: Unit = "days", ) -> float: ... - def two_part_julian_date(self) -> tuple[float, float]: ... - def scale(self) -> Scale: ... - def year(self) -> int: ... - def month(self) -> int: ... - def day(self) -> int: ... - def day_of_year(self) -> int: ... - def hour(self) -> int: ... - def minute(self) -> int: ... - def second(self) -> int: ... - def millisecond(self) -> int: ... - def microsecond(self) -> int: ... - def nanosecond(self) -> int: ... - def picosecond(self) -> int: ... - def femtosecond(self) -> int: ... - def decimal_seconds(self) -> float: ... - def to_tai(self, provider: UT1Provider | None = None) -> Self: ... - def to_tcb(self, provider: UT1Provider | None = None) -> Self: ... - def to_tcg(self, provider: UT1Provider | None = None) -> Self: ... - def to_tdb(self, provider: UT1Provider | None = None) -> Self: ... - def to_tt(self, provider: UT1Provider | None = None) -> Self: ... - def to_ut1(self, provider: UT1Provider | None = None) -> Self: ... - def to_utc(self, provider: UT1Provider | None = None) -> UTC: ... - class TimeDelta: def __new__(cls, seconds: float): ... - def __repr__(self) -> str: ... - def __str__(self) -> str: ... - def __float__(self) -> float: ... - def __neg__(self) -> Self: ... - def __add__(self, other: Self) -> Self: ... - def __sub__(self, other: Self) -> Self: ... - def seconds(self) -> int: ... - def subsecond(self) -> float: ... - @classmethod def from_seconds(cls, seconds: int) -> Self: ... - @classmethod def from_minutes(cls, minutes: float) -> Self: ... - @classmethod def from_hours(cls, hours: float) -> Self: ... - @classmethod def from_days(cls, days: float) -> Self: ... - @classmethod def from_julian_years(cls, years: float) -> Self: ... - @classmethod def from_julian_centuries(cls, centuries: float) -> Self: ... - def to_decimal_seconds(self) -> float: ... - @classmethod def range(cls, start: int, end: int, step: int | None = None) -> list[Self]: ... - class UTC: def __new__( - cls, - year: int, - month: int, - day: int, - hour: int = 0, - minute: int = 0, - seconds: float = 0.0, + cls, + year: int, + month: int, + day: int, + hour: int = 0, + minute: int = 0, + seconds: float = 0.0, ): ... - @classmethod def from_iso(cls, iso: str) -> Self: ... - def __str__(self) -> str: ... - def __repr__(self) -> str: ... - def __eq__(self, other: object) -> bool: ... - def year(self) -> int: ... - def month(self) -> int: ... - def day(self) -> int: ... - def hour(self) -> int: ... - def minute(self) -> int: ... - def second(self) -> int: ... - def millisecond(self) -> int: ... - def microsecond(self) -> int: ... - def nanosecond(self) -> int: ... - def picosecond(self) -> int: ... - def decimal_seconds(self) -> float: ... - def to_tai(self) -> Time: ... - def to_tcb(self) -> Time: ... - def to_tcg(self) -> Time: ... - def to_tdb(self) -> Time: ... - def to_tt(self) -> Time: ... - def to_ut1(self, provider: UT1Provider) -> Time: ... - class UT1Provider: def __new__(cls, path: str): ... - type Origin = Sun | Barycenter | Planet | Satellite | MinorBody - class Sun: def __new__(cls): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def equatorial_radius(self) -> float: ... - class Barycenter: def __new__(cls, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - class Planet: def __new__(cls, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def equatorial_radius(self) -> float: ... - class Satellite: def __new__(cls, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def subplanetary_radius(self) -> float: ... - def along_orbit_radius(self) -> float: ... - class MinorBody: def __new__(cls, name: str): ... - def id(self) -> int: ... - def name(self) -> str: ... - def gravitational_parameter(self) -> float: ... - def mean_radius(self) -> float: ... - def polar_radius(self) -> float: ... - def subplanetary_radius(self) -> float: ... - def along_orbit_radius(self) -> float: ... diff --git a/crates/lox-space/tests/test_states.py b/crates/lox-space/tests/test_states.py index 0a874aab..441f3b15 100644 --- a/crates/lox-space/tests/test_states.py +++ b/crates/lox-space/tests/test_states.py @@ -5,7 +5,18 @@ # file, you can obtain one at https://mozilla.org/MPL/2.0/. import lox_space as lox +import numpy as np +import numpy.testing as npt import pytest +from pathlib import Path + + +@pytest.fixture +def ephemeris(): + spk = ( + Path(__file__).parent.joinpath("..", "..", "..", "data", "de440s.bsp").resolve() + ) + return lox.SPK(str(spk)) def test_state_to_ground_location(): @@ -17,3 +28,31 @@ def test_state_to_ground_location(): assert ground.longitude() == pytest.approx(2.646276127963636) assert ground.latitude() == pytest.approx(-0.2794495715104036) assert ground.altitude() == pytest.approx(417.8524158044338) + + +def test_state_to_origin(ephemeris): + r_venus = np.array( + [ + 1.001977553295792e8, + 2.200234656010247e8, + 9.391473630346918e7, + ] + ) + v_venus = np.array([-59.08617935009049, 22.682387107225292, 12.05029567478702]) + r = np.array([6068279.27, -1692843.94, -2516619.18]) / 1e3 + + v = np.array([-660.415582, 5495.938726, -5303.093233]) / 1e3 + + r_exp = r - r_venus + v_exp = v - v_venus + utc = lox.UTC.from_iso("2016-05-30T12:00:00.000") + tai = utc.to_tai() + + s_earth = lox.State(tai, tuple(r), tuple(v)) + s_venus = s_earth.to_origin(lox.Planet("Venus"), ephemeris) + + r_act = s_venus.position() + v_act = s_venus.velocity() + + npt.assert_allclose(r_act, r_exp) + npt.assert_allclose(v_act, v_exp)