diff --git a/CHANGELOG.md b/CHANGELOG.md index b39eab1..f0e7a95 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - Added worked example of calculating if an object is an Earth Trojan. +- Add Earth obliquity calculation. +- Exposed WGS84 coordinate conversions in the python frontend. ### Fixed diff --git a/src/kete/conversion.py b/src/kete/conversion.py index 03c80ef..36ebd92 100644 --- a/src/kete/conversion.py +++ b/src/kete/conversion.py @@ -9,6 +9,27 @@ from numpy.typing import NDArray from . import _core from . import constants +from ._core import compute_obliquity + +__all__ = [ + "bin_data", + "compute_airmass", + "compute_albedo", + "compute_aphelion", + "compute_diameter", + "compute_earth_radius", + "compute_eccentric_anomaly", + "compute_h_mag", + "compute_obliquity", + "compute_semi_major", + "compute_tisserand", + "dec_degrees_to_dms", + "dec_dms_to_degrees", + "flux_to_mag", + "mag_to_flux", + "ra_degrees_to_hms", + "ra_hms_to_degrees", +] logger = logging.getLogger(__name__) diff --git a/src/kete/rust/frame.rs b/src/kete/rust/frame.rs index 6324f55..4aecd2e 100644 --- a/src/kete/rust/frame.rs +++ b/src/kete/rust/frame.rs @@ -75,3 +75,37 @@ pub fn frame_change_py( pub fn wgs_lat_lon_to_ecef(lat: f64, lon: f64, h: f64) -> (f64, f64, f64) { geodetic_lat_lon_to_ecef(lat.to_radians(), lon.to_radians(), h) } + +/// Compute WCS84 Geodetic latitude/longitude/height from a ECEF position. +/// +/// This returns the lat, lon, and height from the WGS84 oblate Earth. +/// +/// Parameters +/// ---------- +/// x : +/// ECEF x position in km. +/// y : +/// ECEF y position in km. +/// z : +/// ECEF z position in km. +#[pyfunction] +#[pyo3(name = "ecef_to_wgs_lat_lon")] +pub fn ecef_to_wgs_lat_lon(x: f64, y: f64, z: f64) -> (f64, f64, f64) { + ecef_to_geodetic_lat_lon(x, y, z) +} + +/// Calculate the obliquity angle of the Earth at the specified time. +/// +/// This is only valid for several centuries near J2000. +/// +/// The equation is from the 2010 Astronomical Almanac. +/// +/// Parameters +/// ---------- +/// time: +/// Calculate the obliquity angle of the Earth at the specified time. +#[pyfunction] +#[pyo3(name = "compute_obliquity")] +pub fn calc_obliquity_py(time: f64) -> f64 { + calc_obliquity(time).to_degrees() +} diff --git a/src/kete/rust/lib.rs b/src/kete/rust/lib.rs index 6b1fa7a..fe53abf 100644 --- a/src/kete/rust/lib.rs +++ b/src/kete/rust/lib.rs @@ -82,6 +82,8 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(frame::frame_change_py, m)?)?; m.add_function(wrap_pyfunction!(frame::wgs_lat_lon_to_ecef, m)?)?; + m.add_function(wrap_pyfunction!(frame::ecef_to_wgs_lat_lon, m)?)?; + m.add_function(wrap_pyfunction!(frame::calc_obliquity_py, m)?)?; m.add_function(wrap_pyfunction!(kepler::compute_eccentric_anomaly_py, m)?)?; m.add_function(wrap_pyfunction!(kepler::propagation_kepler_py, m)?)?; diff --git a/src/kete/rust/spice/pck.rs b/src/kete/rust/spice/pck.rs index 435f297..0e1178f 100644 --- a/src/kete/rust/spice/pck.rs +++ b/src/kete/rust/spice/pck.rs @@ -1,10 +1,9 @@ -use kete_core::prelude::*; +use kete_core::frames::ecef_to_geodetic_lat_lon; use kete_core::spice::{get_pck_singleton, get_spk_singleton}; +use kete_core::{constants, prelude::*}; use pyo3::{pyfunction, PyResult}; -use crate::frame::PyFrames; use crate::state::PyState; -use crate::vector::Vector; /// Load all specified files into the PCK shared memory singleton. #[pyfunction] @@ -64,29 +63,35 @@ pub fn pck_earth_frame_py( Ok(PyState(state)) } -/// Convert a [`State`] to the Earth's surface reference frame. +/// Convert a [`State`] to the Earth's surface lat/lon/height on the WGS84 reference. /// /// This requires the `earth_000101_*.pck` file to be loaded which contains the /// instantaneous earth frame information. The one provided by kete has dates from -/// around 2000 to early 2024. New files may be downloaded from the NAIF website for -/// additional precision or years of epoch. The current file is accurate to ~5 cm -/// precision. +/// around 2000 to early 2024, along with predicts into the future. New files may be +/// downloaded from the NAIF website for additional precision or years of epoch. /// /// Parameters /// ---------- /// state: State -/// Convert a given state to an geocentered State. +/// Convert the given state to latitude, longitude, and height in km on the WGS84 +/// reference. #[pyfunction] -#[pyo3(name = "pck_state_to_frame")] -pub fn pck_state_to_earth(state: PyState) -> PyResult { +#[pyo3(name = "state_to_earth_pos")] +pub fn pck_state_to_earth(state: PyState) -> PyResult<(f64, f64, f64)> { let pcks = get_pck_singleton().try_read().unwrap(); let state = state.change_center(399)?.as_ecliptic()?; let frame = pcks.try_get_orientation(3000, state.jd())?; let mut state = state.0; state.try_change_frame_mut(frame)?; + let [x, y, z] = state.pos; + let (lat, lon, height) = ecef_to_geodetic_lat_lon( + x * constants::AU_KM, + y * constants::AU_KM, + z * constants::AU_KM, + ); - Ok(Vector::new(state.pos, PyFrames::Undefined)) + Ok((lat.to_degrees(), lon.to_degrees(), height)) } /// Reset the contents of the PCK shared memory to the default set of PCK kernels. diff --git a/src/kete/spice.py b/src/kete/spice.py index 241cc66..95d1bf8 100644 --- a/src/kete/spice.py +++ b/src/kete/spice.py @@ -10,6 +10,7 @@ from .constants import AU_KM from .cache import download_file, cache_path from .vector import Frames, State +from ._core import state_to_earth_pos __all__ = [ "SpkInfo", @@ -23,6 +24,7 @@ "kernel_header_comments", "mpc_code_to_ecliptic", "earth_pos_to_ecliptic", + "state_to_earth_pos", "moon_illumination_frac", ] diff --git a/src/kete/vector.py b/src/kete/vector.py index 767f2b9..a3f3c85 100644 --- a/src/kete/vector.py +++ b/src/kete/vector.py @@ -11,10 +11,20 @@ State, CometElements, SimultaneousStates, + wgs_lat_lon_to_ecef, + ecef_to_wgs_lat_lon, ) -__all__ = ["Frames", "Vector", "State", "CometElements", "SimultaneousStates"] +__all__ = [ + "Frames", + "Vector", + "State", + "CometElements", + "SimultaneousStates", + "wgs_lat_lon_to_ecef", + "ecef_to_wgs_lat_lon", +] Vector.dec_dms = property( fget=lambda self: conversion.dec_degrees_to_dms(self.dec), diff --git a/src/kete_core/src/frames/definitions.rs b/src/kete_core/src/frames/definitions.rs index 6925124..516efc1 100644 --- a/src/kete_core/src/frames/definitions.rs +++ b/src/kete_core/src/frames/definitions.rs @@ -9,6 +9,7 @@ use std::f64::consts::{PI, TAU}; use std::fmt::{self, Debug, Display}; use crate::prelude::{Error, NeosResult}; +use crate::time::Time; /// Coordinate frames. /// @@ -107,6 +108,22 @@ impl Display for Frame { /// - https://ssd.jpl.nasa.gov/horizons/manual.html#defs const OBLIQUITY: f64 = 0.40909280422232897; +/// Compute the angle of obliquity of Earth. +/// +/// This is only valid for several centuries near J2000. +/// +/// The equation here is from the 2010 Astronomical Almanac. +/// +pub fn calc_obliquity(jd: f64) -> f64 { + // centuries from j2000 + let c = (jd - Time::j2000().jd) / 365.25 / 100.0; + (23.439279444444444 + + c * (-0.013010213611111 + + c * (-5.08611111111111e-08 + + c * (5.565e-07 - c * (1.6e-10 + -1.1777777777777779e-11 * c))))) + .to_radians() +} + lazy_static! { static ref ECLIPTIC_EQUATORIAL_ROT: Rotation3 = { let x = nalgebra::Unit::new_unchecked(Vector3::new(1.0, 0.0, 0.0)); diff --git a/src/kete_core/src/spice/daf.rs b/src/kete_core/src/spice/daf.rs index 2463809..b3a4249 100644 --- a/src/kete_core/src/spice/daf.rs +++ b/src/kete_core/src/spice/daf.rs @@ -1,14 +1,14 @@ -/// Support for arbitrary DAF files -/// DAF is a superset which includes SPK and PCK files. -/// -/// DAF files are laid out in 1024 Byte "Records" -/// - The first record is header information about the contents of the file. -/// - The following N records are text comments. -/// - Immediately following the comments there is a Summary Record. -/// -/// These summary records contain the location information for all the contents -/// of the DAF file. -/// +//! Support for arbitrary DAF files +//! DAF is a superset which includes SPK and PCK files. +//! +//! DAF files are laid out in 1024 Byte "Records" +//! - The first record is header information about the contents of the file. +//! - The following N records are text comments. +//! - Immediately following the comments there is a Summary Record. +//! +//! These summary records contain the location information for all the contents +//! of the DAF file. +//! use super::binary::{ bytes_to_f64, bytes_to_f64_vec, bytes_to_i32, bytes_to_i32_vec, bytes_to_string, read_bytes_exact, read_f64_vec, read_str, diff --git a/src/kete_core/src/spice/naif_ids.rs b/src/kete_core/src/spice/naif_ids.rs index a6a7aa1..cc5e7b7 100644 --- a/src/kete_core/src/spice/naif_ids.rs +++ b/src/kete_core/src/spice/naif_ids.rs @@ -1,6 +1,6 @@ -/// List of NAIF ID values. -/// This list is not comprehensive, but is more complete than the C-SPICE -/// implementation. +//! List of NAIF ID values. +//! This list is not comprehensive, but is more complete than the C-SPICE +//! implementation. use lazy_static::lazy_static; use serde::Deserialize; diff --git a/src/kete_core/src/spice/pck.rs b/src/kete_core/src/spice/pck.rs index 9ba3a23..8e9c55a 100644 --- a/src/kete_core/src/spice/pck.rs +++ b/src/kete_core/src/spice/pck.rs @@ -1,12 +1,10 @@ -/// Loading and reading of states from JPL PCK kernel files. -/// -/// PCKs are intended to be loaded into a singleton which is accessible via the -/// [`get_pck_singleton`] function defined below. This singleton is wrapped in a RwLock, -/// meaning before its use it must by unwrapped. A vast majority of intended use cases -/// will only be the read case. -/// ``` -/// -/// +//! Loading and reading of states from JPL PCK kernel files. +//! +//! PCKs are intended to be loaded into a singleton which is accessible via the +//! [`get_pck_singleton`] function defined below. This singleton is wrapped in a RwLock, +//! meaning before its use it must by unwrapped. A vast majority of intended use cases +//! will only be the read case. +//! use super::daf::{DAFType, DafFile}; use super::pck_segments::PckSegment; use crate::errors::{Error, NeosResult}; diff --git a/src/kete_core/src/spice/pck_segments.rs b/src/kete_core/src/spice/pck_segments.rs index 2c800bb..4852dc7 100644 --- a/src/kete_core/src/spice/pck_segments.rs +++ b/src/kete_core/src/spice/pck_segments.rs @@ -1,19 +1,19 @@ +//! Most users should interface with `pck.rs`, not this module. +//! +//! PCK Files are collections of `Segments`, which are ranges of times where the state +//! of an object is recorded. These segments are typically made up of many individual +//! `Records`, with an associated maximum and minimum time where they are valid for. +//! +//! There are unique structs for each possible segment type, not all are currently +//! supported. Each segment type must implement the PCKSegment trait, which allows for +//! the loading and querying of states contained within. +//! +//! +//! +//! There is a lot of repetition in this file, as many of the segment types have very +//! similar internal structures. +//! use super::daf::DafArray; -/// Most users should interface with `pck.rs`, not this module. -/// -/// PCK Files are collections of `Segments`, which are ranges of times where the state -/// of an object is recorded. These segments are typically made up of many individual -/// `Records`, with an associated maximum and minimum time where they are valid for. -/// -/// There are unique structs for each possible segment type, not all are currently -/// supported. Each segment type must implement the PCKSegment trait, which allows for -/// the loading and querying of states contained within. -/// -/// -/// -/// There is a lot of repetition in this file, as many of the segment types have very -/// similar internal structures. -/// use super::interpolation::*; use super::{jd_to_spice_jd, spice_jds_to_jd}; use crate::errors::Error; diff --git a/src/kete_core/src/spice/spk.rs b/src/kete_core/src/spice/spk.rs index c8d1e6b..a331d9d 100644 --- a/src/kete_core/src/spice/spk.rs +++ b/src/kete_core/src/spice/spk.rs @@ -1,23 +1,23 @@ -/// Loading and reading of states from JPL SPK kernel files. -/// -/// SPKs are intended to be loaded into a singleton which is accessible via the -/// [`get_spk_singleton`] function defined below. This singleton is wrapped in a RwLock, -/// meaning before its use it must by unwrapped. A vast majority of intended use cases -/// will only be the read case. -/// -/// Here is a small worked example: -/// ``` -/// use kete_core::spice::get_spk_singleton; -/// use kete_core::frames::Frame; -/// -/// // get a read-only reference to the [`SegmentCollection`] -/// let singleton = get_spk_singleton().try_read().unwrap(); -/// -/// // get the state of 399 (Earth) with respect to the Sun (10) -/// let state = singleton.try_get_state(399, 2451545.0, 10, Frame::Ecliptic); -/// ``` -/// -/// +//! Loading and reading of states from JPL SPK kernel files. +//! +//! SPKs are intended to be loaded into a singleton which is accessible via the +//! [`get_spk_singleton`] function defined below. This singleton is wrapped in a RwLock, +//! meaning before its use it must by unwrapped. A vast majority of intended use cases +//! will only be the read case. +//! +//! Here is a small worked example: +//! ``` +//! use kete_core::spice::get_spk_singleton; +//! use kete_core::frames::Frame; +//! +//! // get a read-only reference to the [`SegmentCollection`] +//! let singleton = get_spk_singleton().try_read().unwrap(); +//! +//! // get the state of 399 (Earth) with respect to the Sun (10) +//! let state = singleton.try_get_state(399, 2451545.0, 10, Frame::Ecliptic); +//! ``` +//! +//! use super::daf::DafFile; use super::{spk_segments::*, DAFType}; use crate::errors::Error; diff --git a/src/kete_core/src/spice/spk_segments.rs b/src/kete_core/src/spice/spk_segments.rs index 703a565..ec9714e 100644 --- a/src/kete_core/src/spice/spk_segments.rs +++ b/src/kete_core/src/spice/spk_segments.rs @@ -1,17 +1,17 @@ -/// Most users should interface with `spk.rs`, not this module. -/// -/// SPK Files are collections of `Segments`, which are ranges of times where the state -/// of an object is recorded. These segments are typically made up of many individual -/// `Records`, with an associated maximum and minimum time where they are valid for. -/// -/// There are unique structs for each possible segment type, not all are currently -/// supported. Each segment type must implement the SPKSegment trait, which allows for -/// the loading and querying of states contained within. -/// -/// -/// -/// There is a lot of repetition in this file, as many of the segment types have very -/// similar internal structures. +//! Most users should interface with `spk.rs`, not this module. +//! +//! SPK Files are collections of `Segments`, which are ranges of times where the state +//! of an object is recorded. These segments are typically made up of many individual +//! `Records`, with an associated maximum and minimum time where they are valid for. +//! +//! There are unique structs for each possible segment type, not all are currently +//! supported. Each segment type must implement the SPKSegment trait, which allows for +//! the loading and querying of states contained within. +//! +//! +//! +//! There is a lot of repetition in this file, as many of the segment types have very +//! similar internal structures. use super::interpolation::*; use super::{jd_to_spice_jd, spice_jds_to_jd, DafArray}; use crate::constants::AU_KM; diff --git a/src/kete_core/src/time/mod.rs b/src/kete_core/src/time/mod.rs index db9448b..d7f1834 100644 --- a/src/kete_core/src/time/mod.rs +++ b/src/kete_core/src/time/mod.rs @@ -158,6 +158,11 @@ impl Time { let datetime = self.to_datetime()?; Ok(datetime.to_rfc3339()) } + + /// J2000 reference time. + pub fn j2000() -> Time { + Time::::new(2451545.0) + } } impl Time {