diff --git a/CHANGELOG.md b/CHANGELOG.md index e082b50..097647d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- Added support for long term integrations, on the scale of less than a mega-year. + ### Changed - Optimizations to SPICE kernel queries leading to a 20% speedup in orbit propagation, @@ -16,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Moved python wrappers over propagation into the rust backend entirely. - Simplified polymorphic support in the rust wrappers, removing some intermediate objects which were over-complicating the interface. +- Refactored the computation of gravitational forces, improving maintainability. ### Fixed diff --git a/docs/conf.py b/docs/conf.py index a531210..d6d987d 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -63,7 +63,7 @@ # -- Sphinx gallery settings -------------------------------------------------- sphinx_gallery_conf = { - "run_stale_examples": False, + "run_stale_examples": True, "filename_pattern": "", "examples_dirs": ["../src/examples"], } diff --git a/src/examples/README.rst b/src/examples/README.rst index 8d9cc91..b399f20 100644 --- a/src/examples/README.rst +++ b/src/examples/README.rst @@ -1,4 +1,4 @@ Examples ======== -Here are a collection of examples which demonstrate different parts of the neospy. +A collection of examples which demonstrate different parts of the neospy. diff --git a/src/neospy/propagation.py b/src/neospy/propagation.py index f340581..6cb30db 100644 --- a/src/neospy/propagation.py +++ b/src/neospy/propagation.py @@ -10,11 +10,17 @@ from .vector import State from . import spice -from ._core import NonGravModel, propagate_n_body, propagate_two_body +from ._core import ( + NonGravModel, + propagate_n_body, + propagate_n_body_long, + propagate_two_body, +) __all__ = [ "propagate_n_body", + "propagate_n_body_long", "propagate_two_body", "NonGravModel", "moid", diff --git a/src/neospy/rust/lib.rs b/src/neospy/rust/lib.rs index 679e63a..d1e9ffc 100644 --- a/src/neospy/rust/lib.rs +++ b/src/neospy/rust/lib.rs @@ -62,6 +62,7 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_function(wrap_pyfunction!(kepler::propagation_kepler_py, m)?)?; m.add_function(wrap_pyfunction!(propagation::propagation_n_body_spk_py, m)?)?; + m.add_function(wrap_pyfunction!(propagation::propagation_n_body_py, m)?)?; m.add_function(wrap_pyfunction!(fovs::fov_checks_py, m)?)?; m.add_function(wrap_pyfunction!(fovs::fov_spk_checks_py, m)?)?; diff --git a/src/neospy/rust/propagation.rs b/src/neospy/rust/propagation.rs index e3dff2b..89da4a7 100644 --- a/src/neospy/rust/propagation.rs +++ b/src/neospy/rust/propagation.rs @@ -1,7 +1,7 @@ use itertools::Itertools; use neospy_core::{ errors::Error, - propagation, + propagation::{self, NonGravModel}, spice::{self, get_spk_singleton}, state::State, time::{scales::TDB, Time}, @@ -136,3 +136,97 @@ pub fn propagation_n_body_spk_py( Ok(res) } + +/// It is *STRONGLY* recommended to use `propagate_n_body` instead of this function +/// wherever possible. This function is specifically meant for kilo-year or longer +/// simulations, it is slower and less accurate than `propagate_n_body`, but that +/// function only works for as long as there are SPICE kernels for the planets +/// available. +/// +/// This is designed to not require SPICE kernels to function, and is meant for long +/// term simulations on the other of less than a mega-year. It is not recommended +/// for orbits longer than this. +/// +/// Propagation using this will treat the Earth and Moon as a single object for +/// performance reasons. +/// +/// Propagate the provided :class:`~neospy.State` using N body mechanics to the +/// specified times, very few approximations are made, this can be very CPU intensive. +/// +/// This does not compute light delay, however it does include corrections for general +/// relativity due to the Sun. +/// +/// This returns two lists of states: +/// - First one contains the states of the objects at the end of the integration +/// - Second contains the states of the planets at the end of the integration, which may +/// be used as input for continuing the integration. +/// +/// Parameters +/// ---------- +/// states: +/// The initial states, this is a list of multiple State objects. +/// jd: +/// A JD to propagate the initial states to. +/// planet_states: +/// Optional list of the planet's states at the same time as the provided states. +/// If this is not provided, the planets positions are loaded from the Spice kernels +/// if that information is available. +/// non_gravs: +/// A list of non-gravitational terms for each object. If provided, then every +/// object must have an associated :class:`~NonGravModel`. +/// batch_size: +/// Number of objects to propagate at once with the planets. This is used to break +/// up the simulation for multi-core support. It additionally has effects on the +/// integrator stepsize which is difficult to predict before running. This can be +/// manually tuned for increased performance, it should have no other effects than +/// performance. +/// +/// Returns +/// ------- +/// Iterable +/// A :class:`~neospy.State` at the new time. +#[pyfunction] +#[pyo3(name = "propagate_n_body_long", signature = (states, jd_final, planet_states=None, non_gravs=None, batch_size=10))] +pub fn propagation_n_body_py( + states: Vec, + jd_final: PyTime, + planet_states: Option>, + non_gravs: Option>>, + batch_size: usize, +) -> PyResult<(Vec, Vec)> { + let states: Vec = states.into_iter().map(|x| x.0).collect(); + let planet_states: Option> = + planet_states.map(|s| s.into_iter().map(|x| x.0).collect()); + + let non_gravs = non_gravs.unwrap_or(vec![None; states.len()]); + let non_gravs: Vec> = + non_gravs.into_iter().map(|y| y.map(|z| z.0)).collect(); + + let jd = jd_final.jd(); + let res = states + .into_iter() + .zip(non_gravs.into_iter()) + .collect_vec() + .par_chunks(batch_size) + .map(|chunk| { + let (chunk_state, chunk_nongrav): (Vec, Vec>) = + chunk.iter().cloned().collect(); + + propagation::propagate_n_body_vec(chunk_state, jd, planet_states.clone(), chunk_nongrav) + .map(|(states, planets)| { + ( + states.into_iter().map(PyState::from).collect::>(), + planets.into_iter().map(PyState::from).collect::>(), + ) + }) + }) + .collect::, _>>()?; + + let mut final_states = Vec::new(); + let mut final_planets = Vec::new(); + for (mut state, planet) in res.into_iter() { + final_states.append(&mut state); + final_planets = planet; + } + Ok((final_states, final_planets)) +} diff --git a/src/neospy/rust/spice/spk.rs b/src/neospy/rust/spice/spk.rs index 88f49ae..3c59a35 100644 --- a/src/neospy/rust/spice/spk.rs +++ b/src/neospy/rust/spice/spk.rs @@ -1,5 +1,5 @@ use neospy_core::spice::{get_spk_singleton, try_name_from_id}; -use pyo3::{pyfunction, PyResult}; +use pyo3::{pyfunction, PyResult, Python}; use crate::frame::PyFrames; use crate::state::PyState; @@ -8,12 +8,13 @@ use crate::time::PyTime; /// Load all specified files into the SPK shared memory singleton. #[pyfunction] #[pyo3(name = "spk_load")] -pub fn spk_load_py(filenames: Vec) -> PyResult<()> { +pub fn spk_load_py(py: Python<'_>, filenames: Vec) -> PyResult<()> { let mut singleton = get_spk_singleton().write().unwrap(); if filenames.len() > 100 { eprintln!("Loading {} spk files...", filenames.len()); } for filename in filenames.iter() { + py.check_signals()?; let load = (*singleton).load_file(filename); if let Err(err) = load { eprintln!("{} failed to load. {}", filename, err); diff --git a/src/neospy_core/benches/propagation.rs b/src/neospy_core/benches/propagation.rs index 9761ed4..a1a3319 100644 --- a/src/neospy_core/benches/propagation.rs +++ b/src/neospy_core/benches/propagation.rs @@ -1,4 +1,6 @@ extern crate criterion; +use std::time::Duration; + use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion}; use lazy_static::lazy_static; use neospy_core::prelude::*; @@ -54,6 +56,16 @@ fn prop_n_body_radau(state: State, dt: f64) { propagation::propagate_n_body_spk(state, jd, false, None).unwrap(); } +fn prop_n_body_vec_radau(mut state: State, dt: f64) { + let spk = get_spk_singleton().read().unwrap(); + spk.try_change_center(&mut state, 10).unwrap(); + state.try_change_frame_mut(Frame::Ecliptic).unwrap(); + let states = vec![state.clone(); 100]; + let non_gravs = vec![None; 100]; + let jd = state.jd + dt; + propagation::propagate_n_body_vec(states, jd, None, non_gravs).unwrap(); +} + fn prop_n_body_radau_par(state: State, dt: f64) { let states: Vec = (0..100).map(|_| state.clone()).collect(); let _tmp: Vec = states @@ -87,7 +99,7 @@ pub fn two_body_numeric(c: &mut Criterion) { _ => panic!(), }; twobody_num_group.bench_with_input(BenchmarkId::new("Single", name), &state, |b, s| { - b.iter(|| prop_2_body_radau(s.clone(), black_box(1000.0))) + b.iter(|| prop_2_body_radau(black_box(s.clone()), black_box(1000.0))) }); } } @@ -106,7 +118,7 @@ pub fn n_body_prop(c: &mut Criterion) { _ => panic!(), }; nbody_group.bench_with_input(BenchmarkId::new("Single", name), &state, |b, s| { - b.iter(|| prop_n_body_radau(s.clone(), black_box(1000.0))) + b.iter(|| prop_n_body_radau(black_box(s.clone()), black_box(1000.0))) }); nbody_group.bench_with_input(BenchmarkId::new("Parallel", name), &state, |b, s| { @@ -114,6 +126,24 @@ pub fn n_body_prop(c: &mut Criterion) { }); } } +pub fn n_body_prop_vec(c: &mut Criterion) { + let mut nbody_group = c.benchmark_group("N-Body-Vec"); + + for state in [ + CIRCULAR.clone(), + ELLIPTICAL.clone(), + PARABOLIC.clone(), + HYPERBOLIC.clone(), + ] { + let name = match &state.desig { + Desig::Name(n) => n, + _ => panic!(), + }; + nbody_group.bench_with_input(BenchmarkId::new("Single", name), &state, |b, s| { + b.iter(|| prop_n_body_vec_radau(black_box(s.clone()), black_box(1000.0))) + }); + } +} pub fn two_body_analytic(c: &mut Criterion) { let mut twobody_group = c.benchmark_group("2-Body-Analytic"); @@ -135,6 +165,6 @@ pub fn two_body_analytic(c: &mut Criterion) { } criterion_group!(name=benches; - config = Criterion::default().with_profiler(PProfProfiler::new(100, Output::Flamegraph(None))); - targets=two_body_analytic, n_body_prop, two_body_numeric); + config = Criterion::default().sample_size(30).measurement_time(Duration::from_secs(15)).with_profiler(PProfProfiler::new(100, Output::Flamegraph(None))); + targets=n_body_prop_vec, two_body_analytic, n_body_prop, two_body_numeric); criterion_main!(benches); diff --git a/src/neospy_core/src/constants/gravity.rs b/src/neospy_core/src/constants/gravity.rs new file mode 100644 index 0000000..ae39b52 --- /dev/null +++ b/src/neospy_core/src/constants/gravity.rs @@ -0,0 +1,268 @@ +use nalgebra::Vector3; + +use crate::frames; + +use super::{AU_KM, C_AU_PER_DAY_INV_SQUARED}; + +/// Standard Gravitational Constants of the Sun +/// AU^3 / (Day^2 * Solar Mass) +pub const GMS: f64 = 0.00029591220828411956; + +/// Gaussian gravitational constant, equivalent to sqrt of GMS. +/// AU^(3/2) per (Day sqrt(Solar Mass)) +pub const GMS_SQRT: f64 = 0.01720209894846; + +/// Sun J2 Parameter +/// +/// This paper below a source, however there are several papers which all put +/// the Sun's J2 at 2.2e-7. +/// +/// "Prospects of Dynamical Determination of General Relativity Parameter β and Solar +/// Quadrupole Moment J2 with Asteroid Radar Astronomy" +/// The Astrophysical Journal, 845:166 (5pp), 2017 August 20 +pub const SUN_J2: f64 = 2.2e-7; + +/// Earth J2 Parameter +/// See "Revisiting Spacetrack Report #3" - Final page of appendix. +pub const EARTH_J2: f64 = 0.00108262998905; + +/// Earth J3 Parameter +pub const EARTH_J3: f64 = -0.00000253215306; + +/// Earth J4 Parameter +pub const EARTH_J4: f64 = -0.00000161098761; + +/// Jupiter J2 Parameter +/// +/// "Measurement of Jupiter’s asymmetric gravity field" +/// https://www.nature.com/articles/nature25776 +/// Nature 555, 220-220, 2018 March 8 +pub const JUPITER_J2: f64 = 0.014696572; + +/// Known massive objects +pub const MASSES: &[GravParams] = &[ + // Sun + GravParams { + naif_id: 10, + mass: GMS, + radius: 0.004654758765894654, + }, + // Mercury Barycenter + GravParams { + naif_id: 1, + mass: 1.66012082548908e-07 * GMS, + radius: 1.63139354095098e-05, + }, + // Venus Barycenter + GravParams { + naif_id: 2, + mass: 2.44783828779694e-06 * GMS, + radius: 4.04537843465442e-05, + }, + // Earth + GravParams { + naif_id: 399, + mass: 3.00348961546514e-06 * GMS, + radius: 4.33036684926559e-05, + }, + // Moon + GravParams { + naif_id: 301, + mass: 3.69430335010988e-08 * GMS, + radius: 1.16138016662292e-05, + }, + // Mars Barycenter + GravParams { + naif_id: 4, + mass: 3.22715608291416e-07 * GMS, + radius: 2.27021279387769e-05, + }, + // Jupiter + GravParams { + naif_id: 5, + mass: 0.0009547919099414246 * GMS, + radius: 0.00047789450254521576, + }, + // Saturn Barycenter + GravParams { + naif_id: 6, + mass: 0.00028588567002459455 * GMS, + radius: 0.0004028666966848747, + }, + // Uranus Barycenter + GravParams { + naif_id: 7, + mass: 4.36624961322212e-05 * GMS, + radius: 0.0001708513622580592, + }, + // Neptune Barycenter + GravParams { + naif_id: 8, + mass: 5.15138377265457e-05 * GMS, + radius: 0.0001655371154958558, + }, + // Ceres + GravParams { + naif_id: 20000001, + mass: 4.7191659919436e-10 * GMS, + radius: 6.276827e-6, + }, + // Pallas + GravParams { + naif_id: 20000002, + mass: 1.0259143964958e-10 * GMS, + radius: 3.415824e-6, + }, + // Vesta + GravParams { + naif_id: 20000004, + mass: 1.3026452498655e-10 * GMS, + radius: 3.512082e-6, + }, + // Hygiea + GravParams { + naif_id: 20000010, + mass: 4.3953391300849e-11 * GMS, + radius: 2.894426e-6, + }, + // Interamnia + GravParams { + naif_id: 20000704, + mass: 1.911e-11 * GMS, + radius: 2.219283e-6, + }, +]; + +/// Earth-Moon Barycenter +pub const EM_BARY: GravParams = GravParams { + naif_id: 3, + mass: (3.00348961546514e-06 + 3.69430335010988e-08) * GMS, + // earth - moon distance + radius: 385_000.0 / AU_KM, +}; + +/// Known planet masses, Sun through Neptune, including the Moon +pub const PLANETS: &[GravParams] = &[ + MASSES[0], MASSES[1], MASSES[2], MASSES[3], MASSES[4], MASSES[5], MASSES[6], MASSES[7], + MASSES[8], MASSES[9], +]; + +/// Known planet masses, Sun through Neptune, Excluding the moon +pub const SIMPLE_PLANETS: &[GravParams] = &[ + MASSES[0], MASSES[1], MASSES[2], EM_BARY, MASSES[5], MASSES[6], MASSES[7], MASSES[8], MASSES[9], +]; + +/// Gravitational model for a basic object +#[derive(Debug, Clone, Copy)] +pub struct GravParams { + /// Associated NAIF id + pub naif_id: i64, + + /// Mass of the object in GMS + pub mass: f64, + + /// Radius of the object in AU. + pub radius: f64, +} + +impl GravParams { + /// Add acceleration to the provided accel vector. + #[inline(always)] + pub fn add_acceleration( + self, + accel: &mut Vector3, + rel_pos: &Vector3, + rel_vel: &Vector3, + ) { + // Basic newtonian gravity + let mass = self.mass; + let r3_inv = rel_pos.norm().powi(-3); + *accel -= &(rel_pos * (r3_inv * mass)); + + // Special cases for different objects + match self.naif_id { + 5 => { + let radius = self.radius; + + // GR correction + apply_gr_correction(accel, rel_pos, rel_vel, &mass); + + // J2 correction + let rel_pos_norm_eclip = frames::equatorial_to_ecliptic(&rel_pos.normalize()); + *accel += frames::ecliptic_to_equatorial(&j2_correction( + rel_pos, + &rel_pos_norm_eclip.z, + &radius, + &JUPITER_J2, + &mass, + )); + } + 10 => { + let radius = self.radius; + + // GR correction + apply_gr_correction(accel, rel_pos, rel_vel, &mass); + + // J2 correction + let rel_pos_norm_eclip = frames::equatorial_to_ecliptic(&rel_pos.normalize()); + *accel += frames::ecliptic_to_equatorial(&j2_correction( + rel_pos, + &rel_pos_norm_eclip.z, + &radius, + &SUN_J2, + &mass, + )); + } + 399 => { + *accel += j2_correction( + rel_pos, + &rel_pos.normalize()[2], + &self.radius, + &EARTH_J2, + &mass, + ) + } + _ => { + // no additional forces beyond newtonian mechanics + } + } + } +} + +/// Calculate the effects of the J2 term +/// +/// Z is the z component of the unit vector. +#[inline(always)] +pub fn j2_correction( + rel_pos: &Vector3, + z: &f64, + radius: &f64, + j2: &f64, + mass: &f64, +) -> Vector3 { + let z_squared = 5.0 * z.powi(2); + let coef = j2 * mass * rel_pos.norm_squared().powi(-5) * 1.5 * radius.powi(2); + Vector3::::new( + -rel_pos[0] * coef * (z_squared - 1.0), + -rel_pos[1] * coef * (z_squared - 1.0), + rel_pos[2] * coef * (z_squared - 3.0), + ) +} + +/// Add the effects of general relativistic motion to an acceleration vector +#[inline(always)] +pub fn apply_gr_correction( + accel: &mut Vector3, + rel_pos: &Vector3, + rel_vel: &Vector3, + mass: &f64, +) { + let r_v = 4.0 * rel_pos.dot(rel_vel); + + let rel_v2: f64 = rel_vel.norm_squared(); + let r = rel_pos.norm(); + + let gr_const: f64 = mass * C_AU_PER_DAY_INV_SQUARED * r.powi(-3); + let c: f64 = 4. * mass / r - rel_v2; + *accel += gr_const * (c * rel_pos + r_v * rel_vel); +} diff --git a/src/neospy_core/src/constants/mod.rs b/src/neospy_core/src/constants/mod.rs index 348e34a..a165a17 100644 --- a/src/neospy_core/src/constants/mod.rs +++ b/src/neospy_core/src/constants/mod.rs @@ -1,10 +1,12 @@ //! # Constants //! Constant values, both universal and observatory specific. //! +mod gravity; mod neos; mod universal; mod wise; +pub use gravity::*; pub use neos::*; pub use universal::*; pub use wise::*; diff --git a/src/neospy_core/src/constants/universal.rs b/src/neospy_core/src/constants/universal.rs index fa88a52..98dd2a5 100644 --- a/src/neospy_core/src/constants/universal.rs +++ b/src/neospy_core/src/constants/universal.rs @@ -1,46 +1,9 @@ -//! # Universal Constants - -/// Standard Gravitational Constants of the Sun -/// AU^3 / (Day^2 * Solar Mass) -pub const GMS: f64 = 0.00029591220828411956; - -/// Gaussian gravitational constant, equivalent to sqrt of GMS. -/// AU^(3/2) per (Day sqrt(Solar Mass)) -pub const GMS_SQRT: f64 = 0.01720209894846; - /// effective surface temperature of the sun in K. pub const SUN_TEMP: f64 = 5778.0; /// Sun diameter in km pub const SUN_DIAMETER: f64 = 1.3914e6; -/// Sun J2 Parameter -/// -/// This paper below a source, however there are several papers which all put -/// the Sun's J2 at 2.2e-7. -/// -/// "Prospects of Dynamical Determination of General Relativity Parameter β and Solar -/// Quadrupole Moment J2 with Asteroid Radar Astronomy" -/// The Astrophysical Journal, 845:166 (5pp), 2017 August 20 -pub const SUN_J2: f64 = 2.2e-7; - -/// Earth J2 Parameter -/// See "Revisiting Spacetrack Report #3" - Final page of appendix. -pub const EARTH_J2: f64 = 0.00108262998905; - -/// Earth J3 Parameter -pub const EARTH_J3: f64 = -0.00000253215306; - -/// Earth J4 Parameter -pub const EARTH_J4: f64 = -0.00000161098761; - -/// Jupiter J2 Parameter -/// -/// "Measurement of Jupiter’s asymmetric gravity field" -/// https://www.nature.com/articles/nature25776 -/// Nature 555, 220-220, 2018 March 8 -pub const JUPITER_J2: f64 = 0.014696572; - /// W/m^2 measured at 1 AU. pub const SOLAR_FLUX: f64 = 1360.8; @@ -69,43 +32,5 @@ pub const GOLDEN_RATIO: f64 = 1.618033988749894; /// pub const V_MAG_ZERO: f64 = 3597.28; -/// The ID value for the massive objects which needs to be numerically evaluated for -/// the effects of gravity. -/// -/// Recorded values are: -/// (ID, GM of the object in solar masses, Radius of the object in AU) -pub const MASSIVE_OBJECTS: &[(i64, f64, f64)] = &[ - (10, 1.0, 0.004654758765894654), // Sun - (1, 1.66012082548908e-07, 1.63139354095098e-05), // Mercury Barycenter - (2, 2.44783828779694e-06, 4.04537843465442e-05), // Venus Barycenter - (399, 3.00348961546514e-06, 4.33036684926559e-05), // Earth, not 3 because earth != barycenter - (301, 3.69430335010988e-08, 1.16138016662292e-05), // Moon - (4, 3.2271560829139e-07, 2.27021279387769e-05), // Mars Barycenter - (5, 0.0009547919099414246, 0.00047789450254521576), // Jupiter Barycenter - (6, 0.00028588567002459455, 0.0004028666966848747), // Saturn Barycenter - (7, 4.36624961322212e-05, 0.0001708513622580592), // Uranus Barycenter - (8, 5.15138377265457e-05, 0.0001655371154958558), // Neptune Barycenter -]; - -/// Recorded values are: -/// (ID, GM of the object in solar masses, Radius of the object in AU) -pub const MASSIVE_OBJECTS_EXTENDED: &[(i64, f64, f64)] = &[ - (10, 1.0, 0.004654758765894654), // Sun - (1, 1.66012082548908e-07, 1.63139354095098e-05), // Mercury Barycenter - (2, 2.44783828779694e-06, 4.04537843465442e-05), // Venus Barycenter - (399, 3.00348961546514e-06, 4.33036684926559e-05), // Earth, not 3 because earth != barycenter - (301, 3.69430335010988e-08, 1.16138016662292e-05), // Moon - (4, 3.2271560829139e-07, 2.27021279387769e-05), // Mars Barycenter - (5, 0.0009547919099414246, 0.00047789450254521576), // Jupiter Barycenter - (6, 0.00028588567002459455, 0.0004028666966848747), // Saturn Barycenter - (7, 4.36624961322212e-05, 0.0001708513622580592), // Uranus Barycenter - (8, 5.15138377265457e-05, 0.0001655371154958558), // Neptune Barycenter - (20000001, 4.7191659919436e-10, 6.276827e-6), // Ceres - (20000002, 1.0259143964958e-10, 3.415824e-6), // Pallas - (20000004, 1.3026452498655e-10, 3.512082e-6), // Vesta - (20000010, 4.3953391300849e-11, 2.894426e-6), // Hygiea - (20000704, 1.911e-11, 2.219283e-6), // Interamnia -]; - /// V-band constant for the relationship between D, H_V, and p_v, in km. pub const C_V: f64 = 1329.0; diff --git a/src/neospy_core/src/propagation/acceleration.rs b/src/neospy_core/src/propagation/acceleration.rs index 09c1d44..178ed81 100644 --- a/src/neospy_core/src/propagation/acceleration.rs +++ b/src/neospy_core/src/propagation/acceleration.rs @@ -11,21 +11,18 @@ //! //! Where `x` and its derivative `x_der` are vectors. This also accepts a mutable //! reference to a metadata collection. Metadata may include things like object -//! specific orbit parameters such as the non-grav A terms, or keep track of close -//! +//! specific orbit parameters such as the non-grav terms, or keep track of close //! //! `exact_eval` is a bool which is passed when the integrator is passing values where //! the `x` and `x_der` are being evaluated at true locations. IE: where the integrator //! thinks that the object should actually be. These times are when close encounter //! information should be recorded. //! -use crate::frames; use crate::prelude::NeosResult; use crate::spice::get_spk_singleton; use crate::{constants::*, errors::Error, frames::Frame, propagation::nongrav::NonGravModel}; -use itertools::Itertools; use nalgebra::allocator::Allocator; -use nalgebra::{DefaultAllocator, Dim, Matrix, Matrix3, OMatrix, OVector, Vector3, U1, U2}; +use nalgebra::{DefaultAllocator, Dim, Matrix3, OVector, Vector3, U1, U2}; use std::ops::AddAssign; /// Metadata object used by the [`central_accel`] function below. @@ -89,17 +86,16 @@ pub struct AccelSPKMeta<'a> { /// This records the ID of the object, time, and distance in AU. pub close_approach: Option<(i64, f64, f64)>, - /// `A` terms of the non-gravitational forces. + /// The non-gravitational forces. /// If this is not provided, only standard gravitational model is applied. - /// If these values are provided, then the effects of the A terms are added in - /// addition to standard forces. - pub non_grav_a: Option, + /// If this are provided, then the effects of the Non-Grav terms are added. + pub non_grav_model: Option, /// The list of massive objects to apply during SPK computation. /// This list contains the ID of the object in the SPK along with the mass and /// radius of the object. Mass is given in fractions of solar mass and radius is /// in AU. - pub massive_obj: &'a [(i64, f64, f64)], + pub massive_obj: &'a [GravParams], } /// Compute the accel on an object which experiences acceleration due to all massive @@ -141,84 +137,31 @@ pub fn spk_accel( let spk = get_spk_singleton().try_read().unwrap(); - for (id, mass, radius) in meta.massive_obj.iter() { - let state = spk.try_get_state(*id, time, 0, Frame::Equatorial)?; + for grav_params in meta.massive_obj.iter() { + let id = grav_params.naif_id; + let radius = grav_params.radius; + let state = spk.try_get_state(id, time, 0, Frame::Equatorial)?; let rel_pos: Vector3 = pos - Vector3::from(state.pos); - let rel_pos_norm = rel_pos.normalize(); - let r = rel_pos.norm(); + let rel_vel: Vector3 = vel - Vector3::from(state.vel); if exact_eval { + let r = rel_pos.norm(); if let Some(close_approach) = meta.close_approach.as_mut() { if close_approach.2 == r { - *close_approach = (*id, r, time); + *close_approach = (id, r, time); } } - } - if r <= *radius { - Err(Error::Impact(*id, time))?; - } - - let r2_inv = r.powi(-2); - accel -= rel_pos_norm * (r2_inv * *mass * GMS); - - // Corrections for the Sun - if *id == 10 { - let r3_inv = r.powi(-3); - let rel_vel: Vector3 = vel - Vector3::from(state.vel); - - let r_v = 4.0 * rel_pos.dot(&rel_vel); - - let rel_v2: f64 = rel_vel.norm_squared(); - - // GR Correction - let gr_const: f64 = mass * C_AU_PER_DAY_INV_SQUARED * r3_inv * GMS; - let c: f64 = 4. * mass * GMS / r - rel_v2; - accel += gr_const * (c * rel_pos + r_v * rel_vel); - - // J2 for the Sun - let coef = SUN_J2 * GMS * *mass * r.powi(-5) * 1.5 * radius.powi(2); - let rel_pos_norm_eclip = frames::equatorial_to_ecliptic(&rel_pos_norm); - let z2 = 5.0 * rel_pos_norm_eclip.z.powi(2); - accel[0] -= rel_pos.x * coef * (z2 - 1.0); - accel[1] -= rel_pos.y * coef * (z2 - 1.0); - accel[2] -= rel_pos.z * coef * (z2 - 3.0); - // non-gravitational forces - if let Some(model) = &meta.non_grav_a { - accel += model.accel_vec(&rel_pos, &rel_vel) + if r <= radius { + Err(Error::Impact(id, time))?; } + } + grav_params.add_acceleration(&mut accel, &rel_pos, &rel_vel); - // Corrections for Jupiter - } else if *id == 5 { - // GR Correction - let r3_inv = r.powi(-3); - let rel_vel: Vector3 = vel - Vector3::from(state.vel); - - let r_v = 4.0 * rel_pos.dot(&rel_vel); - - let rel_v2: f64 = rel_vel.norm_squared(); - let gr_const: f64 = mass * C_AU_PER_DAY_INV_SQUARED * r3_inv * GMS; - let c: f64 = 4. * mass * GMS / r - rel_v2; - accel += gr_const * (c * rel_pos + r_v * rel_vel); - - // J2 for Jupiter - // https://www.nature.com/articles/nature25776 - // Jupiter's pole is 2.3 degrees off of the ecliptic, below is an approximation of this. - let coef = JUPITER_J2 * GMS * *mass * r.powi(-5) * 1.5 * radius.powi(2); - let rel_pos_norm_eclip = frames::equatorial_to_ecliptic(&rel_pos_norm); - let z2 = 5.0 * rel_pos_norm_eclip.z.powi(2); - accel[0] -= rel_pos.x * coef * (z2 - 1.0); - accel[1] -= rel_pos.y * coef * (z2 - 1.0); - accel[2] -= rel_pos.z * coef * (z2 - 3.0); - - // Corrections for Earth - } else if *id == 399 { - // J2 for the Earth - let coef = EARTH_J2 * GMS * *mass * r.powi(-5) * 1.5 * radius.powi(2); - let z2 = 5.0 * rel_pos_norm.z.powi(2); - accel[0] -= rel_pos.x * coef * (z2 - 1.0); - accel[1] -= rel_pos.y * coef * (z2 - 1.0); - accel[2] -= rel_pos.z * coef * (z2 - 3.0); + if grav_params.naif_id == 10 { + if let Some(non_grav) = &meta.non_grav_model { + non_grav.add_acceleration(&mut accel, &rel_pos, &rel_vel) + } } } Ok(accel) @@ -226,21 +169,17 @@ pub fn spk_accel( /// Metadata for the [`vec_accel`] function defined below. #[derive(Debug)] -pub struct AccelVecMeta<'a, D: Dim> -where - DefaultAllocator: Allocator, -{ - /// `A` terms of the non-gravitational forces. +pub struct AccelVecMeta<'a> { + /// The non-gravitational forces. /// If this is not provided, only standard gravitational model is applied. - /// If these values are provided, then the effects of the A terms are added in - /// addition to standard forces. - pub non_grav_a: Option>, + /// If these values are provided, then the effects of the Non-Grav terms are added. + pub non_gravs: Vec>, /// The list of massive objects to apply during SPK computation. /// This list contains the ID of the object in the SPK along with the mass and /// radius of the object. Mass is given in fractions of solar mass and radius is /// in AU. - pub massive_obj: &'a [(i64, f64, f64)], + pub massive_obj: &'a [GravParams], } /// Compute the accel on an object which experiences acceleration due to all massive @@ -259,8 +198,8 @@ pub fn vec_accel( time: f64, pos: &OVector, vel: &OVector, - meta: &mut AccelVecMeta, - _exact_eval: bool, + meta: &mut AccelVecMeta, + exact_eval: bool, ) -> NeosResult> where DefaultAllocator: Allocator + Allocator, @@ -269,64 +208,43 @@ where // (x, y, z, x, y, z, ... let n_objects = pos.len() / 3; + let n_massive = meta.massive_obj.len(); let (dim, _) = pos.shape_generic(); - let mut accel = Matrix::zeros_generic(dim, U1); + let mut accel = OVector::::zeros_generic(dim, U1); + + let mut accel_working = Vector3::zeros(); for idx in 0..n_objects { - let pos_idx = pos.rows(idx * 3, 3); + let pos_idx = pos.fixed_rows::<3>(idx * 3); + let vel_idx = vel.fixed_rows::<3>(idx * 3); - for (idy, (id, mass, radius)) in meta.massive_obj.iter().enumerate() { + for (idy, grav_params) in meta.massive_obj.iter().enumerate() { if idx == idy { continue; } - let pos_idy = pos.rows(idy * 3, 3); + accel_working.fill(0.0); + let radius = grav_params.radius; + let pos_idy = pos.fixed_rows::<3>(idy * 3); + let vel_idy = vel.fixed_rows::<3>(idy * 3); + let rel_pos = pos_idx - pos_idy; - let rel_pos_norm = rel_pos.normalize(); - let r = rel_pos.norm(); + let rel_vel = vel_idx - vel_idy; - if r <= *radius { - Err(Error::Impact(*id, time))?; + if exact_eval & (rel_pos.norm() <= radius) { + Err(Error::Impact(grav_params.naif_id, time))?; } + grav_params.add_acceleration(&mut accel_working, &rel_pos, &rel_vel); - let r2_inv = r.powi(-2); - accel - .rows_mut(idx * 3, 3) - .add_assign(-&rel_pos_norm * r2_inv * *mass * GMS); - - // if it is the sun or jupiter, apply a correction for GR - if *id == 10 || *id == 5 { - let vel_idx = vel.rows(idx * 3, 3); - let vel_idy = vel.rows(idy * 3, 3); - let r3_inv = r.powi(-3); - let rel_vel = vel_idx - vel_idy; - - let r_v = 4.0 * rel_pos.dot(&rel_vel); - - let rel_v2: f64 = rel_vel.norm_squared(); - let gr_const: f64 = mass * C_AU_PER_DAY_INV_SQUARED * r3_inv * GMS; - let c: f64 = 4. * mass * GMS / r - rel_v2; - - accel - .rows_mut(idx * 3, 3) - .add_assign(gr_const * (c * rel_pos + r_v * &rel_vel)); - - // Add non-grav forces if defined. - if *id == 10 { - if let Some(a_matrix) = meta.non_grav_a.as_ref() { - let a_row = a_matrix.row(idx); - let (a1, a2) = a_row.iter().collect_tuple().unwrap(); - accel - .rows_mut(idx * 3, 3) - .add_assign(&rel_pos_norm * r2_inv * *mass * GMS * *a1); - accel - .rows_mut(idx * 3, 3) - .add_assign(rel_vel.normalize() * *mass * GMS * *a2); - } + if (grav_params.naif_id == 10) & (idx > n_massive) { + if let Some(non_grav) = &meta.non_gravs[idx - n_massive] { + non_grav.add_acceleration(&mut accel_working, &rel_pos, &rel_vel); } } + accel.fixed_rows_mut::<3>(idx * 3).add_assign(accel_working); } } + Ok(accel) } @@ -384,8 +302,10 @@ mod tests { let mut pos: Vec = Vec::new(); let mut vel: Vec = Vec::new(); - for (id, _mass, _radius) in MASSIVE_OBJECTS.iter() { - let planet = spk.try_get_state(*id, jd, 0, Frame::Equatorial).unwrap(); + for obj in MASSES.iter() { + let planet = spk + .try_get_state(obj.naif_id, jd, 0, Frame::Equatorial) + .unwrap(); pos.append(&mut planet.pos.into()); vel.append(&mut planet.vel.into()); } @@ -398,15 +318,15 @@ mod tests { &pos.into(), &vel.into(), &mut AccelVecMeta { - non_grav_a: None, - massive_obj: MASSIVE_OBJECTS, + non_gravs: vec![None], + massive_obj: MASSES, }, false, ) .unwrap() .iter() .copied() - .skip(MASSIVE_OBJECTS.len() * 3) + .skip(MASSES.len() * 3) .collect_vec(); let accel2 = spk_accel( @@ -415,8 +335,8 @@ mod tests { &[0.0, 0.0, 1.0].into(), &mut AccelSPKMeta { close_approach: None, - non_grav_a: None, - massive_obj: MASSIVE_OBJECTS, + non_grav_model: None, + massive_obj: MASSES, }, false, ) diff --git a/src/neospy_core/src/propagation/mod.rs b/src/neospy_core/src/propagation/mod.rs index 2cf5130..984f872 100644 --- a/src/neospy_core/src/propagation/mod.rs +++ b/src/neospy_core/src/propagation/mod.rs @@ -3,7 +3,7 @@ //! There are multiple levels of precision available, each with different pros/cons //! (usually performance related). -use crate::constants::{MASSIVE_OBJECTS, MASSIVE_OBJECTS_EXTENDED}; +use crate::constants::{MASSES, PLANETS, SIMPLE_PLANETS}; use crate::errors::Error; use crate::frames::Frame; use crate::prelude::{Desig, NeosResult}; @@ -71,7 +71,7 @@ pub fn propagate_n_body_spk( mut state: State, jd_final: f64, include_extended: bool, - a_terms: Option, + non_grav_model: Option, ) -> NeosResult { let center = state.center_id; let frame = state.frame; @@ -81,15 +81,15 @@ pub fn propagate_n_body_spk( let mass_list = { if include_extended { - MASSIVE_OBJECTS_EXTENDED + MASSES } else { - MASSIVE_OBJECTS + PLANETS } }; let metadata = AccelSPKMeta { close_approach: None, - non_grav_a: a_terms, + non_grav_model, massive_obj: mass_list, }; @@ -139,42 +139,78 @@ pub fn propagation_central(state: &State, jd_final: f64) -> NeosResult<[[f64; 3] /// Propagate using n-body mechanics but skipping SPK queries. /// This will propagate all planets and the Moon, so it may vary from SPK states slightly. -pub fn propagation_n_body_vec(states: Vec, jd_final: f64) -> NeosResult> { - let spk = get_spk_singleton().try_read().unwrap(); - +pub fn propagate_n_body_vec( + states: Vec, + jd_final: f64, + planet_states: Option>, + non_gravs: Vec>, +) -> NeosResult<(Vec, Vec)> { if states.is_empty() { Err(Error::ValueError( "State vector is empty, propagation cannot continue".into(), ))?; } + + if non_gravs.len() != states.len() { + Err(Error::ValueError( + "Number of non-grav models doesnt match the number of provided objects.".into(), + ))?; + } + let jd_init = states.first().unwrap().jd; let mut pos: Vec = Vec::new(); let mut vel: Vec = Vec::new(); let mut desigs: Vec = Vec::new(); - for (id, _mass, _radius) in MASSIVE_OBJECTS.iter() { - let planet = spk.try_get_state(*id, jd_init, 0, Frame::Ecliptic)?; - pos.append(&mut planet.pos.into()); - vel.append(&mut planet.vel.into()); - desigs.push(planet.desig.to_owned()); + let planet_states = planet_states.unwrap_or_else(|| { + let spk = get_spk_singleton().try_read().unwrap(); + let mut planet_states = Vec::with_capacity(SIMPLE_PLANETS.len()); + for obj in SIMPLE_PLANETS.iter() { + let planet = spk + .try_get_state(obj.naif_id, jd_init, 10, Frame::Equatorial) + .expect("Failed to find state for the provided initial jd"); + planet_states.push(planet); + } + planet_states + }); + + if planet_states.len() != SIMPLE_PLANETS.len() { + Err(Error::ValueError( + "Input planet states must contain the correct number of states.".into(), + ))?; + } + if planet_states.first().unwrap().jd != jd_init { + Err(Error::ValueError( + "Planet states JD must match JD of input state.".into(), + ))?; + } + for planet_state in planet_states.into_iter() { + pos.append(&mut planet_state.pos.into()); + vel.append(&mut planet_state.vel.into()); + desigs.push(planet_state.desig); } for mut state in states.into_iter() { - spk.try_change_center(&mut state, 0)?; - state.try_change_frame_mut(Frame::Ecliptic)?; + state.try_change_frame_mut(Frame::Equatorial)?; if jd_init != state.jd { Err(Error::ValueError( "All input states must have the same JD".into(), ))?; } + if state.center_id != 10 { + Err(Error::ValueError( + "Center of all states must be 10 (the Sun).".into(), + ))? + } pos.append(&mut state.pos.into()); vel.append(&mut state.vel.into()); - desigs.push(state.desig.to_owned()); + desigs.push(state.desig); } + let meta = AccelVecMeta { - non_grav_a: None, - massive_obj: MASSIVE_OBJECTS, + non_gravs, + massive_obj: SIMPLE_PLANETS, }; let (pos, vel, _) = { @@ -187,21 +223,15 @@ pub fn propagation_n_body_vec(states: Vec, jd_final: f64) -> NeosResult = Vec::new(); + let sun_pos = pos.fixed_rows::<3>(0); + let sun_vel = vel.fixed_rows::<3>(0); + let mut all_states: Vec = Vec::new(); for (idx, desig) in desigs.into_iter().enumerate() { - let pos_chunk = pos.rows(idx * 3, 3) - sun_pos; - let vel_chunk = vel.rows(idx * 3, 3) - sun_vel; - let state = State::new( - desig, - jd_final, - nalgebra::convert(pos_chunk), - nalgebra::convert(vel_chunk), - Frame::Ecliptic, - 10, - ); - final_states.push([state.pos, state.vel]); + let pos = pos.fixed_rows::<3>(idx * 3) - sun_pos; + let vel = vel.fixed_rows::<3>(idx * 3) - sun_vel; + let state = State::new(desig, jd_final, pos, vel, Frame::Equatorial, 10); + all_states.push(state); } - Ok(final_states) + let final_states = all_states.split_off(SIMPLE_PLANETS.len()); + Ok((final_states, all_states)) } diff --git a/src/neospy_core/src/propagation/nongrav.rs b/src/neospy_core/src/propagation/nongrav.rs index ca2cf6a..0d403b3 100644 --- a/src/neospy_core/src/propagation/nongrav.rs +++ b/src/neospy_core/src/propagation/nongrav.rs @@ -81,7 +81,7 @@ impl NonGravModel { } /// Construct a new non-grav dust model. - pub fn new_dust_raw(beta: f64) -> Self { + pub fn new_dust(beta: f64) -> Self { Self::Dust { beta } } @@ -101,14 +101,20 @@ impl NonGravModel { /// Compute the non-gravitational acceleration vector when provided the position /// and velocity vector with respect to the sun. - pub fn accel_vec(&self, pos: &Vector3, vel: &Vector3) -> Vector3 { + #[inline(always)] + pub fn add_acceleration( + &self, + accel: &mut Vector3, + pos: &Vector3, + vel: &Vector3, + ) { match self { Self::Dust { beta } => { let pos_norm = pos.normalize(); let r_dot = &pos_norm.dot(vel); let norm2_inv = pos.norm_squared().recip(); let scaling = GMS * beta * norm2_inv; - scaling + *accel += scaling * ((1.0 - r_dot * C_AU_PER_DAY_INV_SQUARED) * pos_norm - vel * C_AU_PER_DAY_INV_SQUARED) } @@ -128,10 +134,9 @@ impl NonGravModel { let scale = alpha * rr0.powf(-m) * (1.0 + rr0.powf(*n)).powf(-k); let t_vec = (vel - pos_norm * vel.dot(&pos_norm)).normalize(); let n_vec = t_vec.cross(&pos_norm).normalize(); - let mut accel = pos_norm * (scale * a1); - accel += t_vec * (scale * a2); - accel += n_vec * (scale * a3); - accel + *accel += pos_norm * (scale * a1); + *accel += t_vec * (scale * a2); + *accel += n_vec * (scale * a3); } } }