Release v1.0.2 (#133)

* Add covariance matrix sampling for the HorizonsProperties

* Add support for time-delayed non-gravitational forces (#132)

* Add python 3.13 to the release

* Release v1.0.2

.github/workflows/release-wheels.yml

@@ -23,7 +23,7 @@ jobs:
       # Ensure that a wheel builder finishes even if another fails
       fail-fast: false
       matrix:
-        python: [39, 310, 311, 312]
+        python: [39, 310, 311, 312, 313]
         os: [windows-latest, ubuntu-latest, macos-latest]
         include:
           # Window 64 bit

CHANGELOG.md

@@ -6,6 +6,16 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
+## [v1.0.2]
+
+### Added
+
+- Added Python 3.13 to the built packages.
+- Added `sample` to the `HorizonsProperties` object, allowing sampling of the orbit's
+  uncertainty.
+- Added support for time delayed non-gravitational forces, as is found a number of
+  comets in JPL Horizons.
+
 ## [v1.0.1]
 
 ### Added

Cargo.toml

@@ -1,7 +1,7 @@
 
 [package]
 name = "_core"
-version = "1.0.1"
+version = "1.0.2"
 edition = "2021"
 readme = "README.md"
 license = "BSD-3-Clause"

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "kete"
-version = "1.0.1"
+version = "1.0.2"
 description = "Kete Asteroid Survey Tools"
 readme = "README.md"
 authors = [{name = "Dar Dahlen", email = "[email protected]"},

src/examples/plot_uncertainty.py

@@ -33,18 +33,7 @@
 jds = np.arange(jd_s, jd_e, time_step)
 
 # Sample the covariance matrix
-cov = np.array(obj.covariance.cov_matrix)
-samples = kete.covariance.generate_sample_from_cov(n_samples, cov)
-labels, values = zip(*obj.covariance.params)
-values = np.array(values)
-
-# Construct new states from the covariance samples.
-# Keep the current best estimate as the first state.
-states = [obj.state]
-for sample in samples:
-    sample = sample + values
-    params = dict(zip(labels, sample))
-    states.append(kete.CometElements(obj.desig, obj.covariance.epoch, **params).state)
+states, _ = obj.sample(n_samples)
 
 # Propagate the position of all states to all time steps, recording the V mags
 mags = []

src/kete/horizons.py

@@ -6,14 +6,34 @@
 import json
 import numpy as np
 
-from ._core import HorizonsProperties, Covariance, NonGravModel
+from ._core import HorizonsProperties, Covariance, NonGravModel, CometElements
 from .mpc import unpack_designation, pack_designation
 from .time import Time
 from .cache import cache_path
+from .covariance import generate_sample_from_cov
 
 __all__ = ["HorizonsProperties", "fetch_spice_kernel"]
 
 
+_PARAM_MAP = {
+    "a1": "a1",
+    "a2": "a2",
+    "a3": "a3",
+    "aln": "alpha",
+    "nm": "m",
+    "r0": "r_0",
+    "nk": "k",
+    "nn": "n",
+    "dt": "dt",
+    "e": "eccentricity",
+    "q": "peri_dist",
+    "tp": "peri_time",
+    "node": "lon_of_ascending",
+    "peri": "peri_arg",
+    "i": "inclination",
+}
+
+
 def fetch(name, update_name=True, cache=True, update_cache=False, exact_name=False):
     """
     Fetch object properties from JPL Horizons service.
@@ -113,6 +133,9 @@ def fetch(name, update_name=True, cache=True, update_cache=False, exact_name=Fal
                 elements = {
                     lab: phys[lookup_rev[lab]] for lab in labels if lab in lookup_rev
                 }
+            if "model_pars" in props["orbit"]:
+                for param in props["orbit"]["model_pars"]:
+                    elements[param["name"]] = float(param["value"])
             params = [(lookup_rev.get(x, x), elements.get(x, np.nan)) for x in labels]
             phys["covariance"] = Covariance(name, cov_epoch, params, mat)
     else:
@@ -228,12 +251,7 @@ def _json(self) -> dict:
     return _fetch_json(self.desig, update_cache=False)
 
 
-@property  # type: ignore
-def _nongrav(self) -> NonGravModel:
-    """
-    The non-gravitational forces model from the values returned from horizons.
-    """
-
+def _nongrav_params(self) -> dict:
     # default parameters used by jpl horizons for their non-grav models.
     params = {
         "a1": 0.0,
@@ -244,34 +262,80 @@ def _nongrav(self) -> NonGravModel:
         "m": 2.15,
         "n": 5.093,
         "k": 4.6142,
-    }
-
-    lookup = {
-        "A1": "a1",
-        "A2": "a2",
-        "A3": "a3",
-        "ALN": "alpha",
-        "NM": "m",
-        "R0": "r_0",
-        "NK": "k",
-        "NN": "n",
+        "dt": 0.0,
     }
 
     orbit = self.json["orbit"]
     if "model_pars" not in orbit:
-        return None
+        return params
 
     for vals in orbit["model_pars"]:
-        if vals["name"] not in lookup:
+        if vals["name"].lower() not in params:
             raise ValueError("Unknown non-grav values: ", vals)
-        params[lookup[vals["name"]]] = float(vals["value"])
+        params[_PARAM_MAP[vals["name"].lower()]] = float(vals["value"])
+    return params
 
+
+@property  # type: ignore
+def _nongrav(self) -> NonGravModel:
+    """
+    The non-gravitational forces model from the values returned from horizons.
+    """
+    params = _nongrav_params(self)
     return NonGravModel.new_comet(**params)
 
 
+def _sample(self, n_samples):
+    """
+    Sample the covariance matrix for this object, returning `n_samples` of states and
+    non-gravitational models, which may be used to propagate the orbit in time.
+
+    This uses the full covariance matrix in order to correctly sample the object's orbit
+    with all correlations, including correlations with the non-gravitational forces.
+
+    Parameters
+    ----------
+    n_samples :
+        The number of samples to take of the covariance.
+    """
+    matrix = self.covariance.cov_matrix
+    epoch = self.covariance.epoch
+    samples = generate_sample_from_cov(n_samples, matrix)
+
+    elem_keywords = [
+        "eccentricity",
+        "inclination",
+        "lon_of_ascending",
+        "peri_arg",
+        "peri_dist",
+        "peri_time",
+    ]
+
+    labels = self.json["orbit"]["covariance"]["labels"]
+    mapped_label = [_PARAM_MAP[k.lower()] for k in labels]
+
+    best_params = _nongrav_params(self)
+    for prop in elem_keywords:
+        best_params[prop] = getattr(self, prop)
+    states = []
+    non_gravs = []
+    for sample in samples:
+        cur_params = {
+            label: best_params[label] + d for label, d in zip(mapped_label, sample)
+        }
+        elem_params = {x: cur_params.pop(x) for x in elem_keywords}
+        state = CometElements(self.desig, epoch, **elem_params).state
+        non_grav = NonGravModel.new_comet(**cur_params)
+        states.append(state)
+        non_gravs.append(non_grav)
+
+    return states, non_gravs
+
+
 HorizonsProperties.fetch = fetch
 HorizonsProperties.json = _json
 HorizonsProperties.non_grav = _nongrav
+HorizonsProperties.sample = _sample
 
 
 def fetch_spice_kernel(

src/kete/rust/nongrav.rs

@@ -97,8 +97,11 @@ impl PyNonGravModel {
     ///
     /// When alpha=1.0, n=0.0, k=0.0, r0=1.0, and m=2.0, this is equivalent to a
     /// :math:`1/r^2` correction.
+    ///
+    /// This includes an optional time delay, which the non-gravitational forces are
+    /// time delayed.
     #[allow(clippy::too_many_arguments)]
-    #[pyo3(signature = (a1, a2, a3, alpha=0.1112620426, r_0=2.808, m=2.15, n=5.093, k=4.6142))]
+    #[pyo3(signature = (a1=0.0, a2=0.0, a3=0.0, alpha=0.1112620426, r_0=2.808, m=2.15, n=5.093, k=4.6142, dt=0.0))]
     #[staticmethod]
     pub fn new_comet(
         a1: f64,
@@ -109,6 +112,7 @@ impl PyNonGravModel {
         m: f64,
         n: f64,
         k: f64,
+        dt: f64,
     ) -> Self {
         Self(NonGravModel::JplComet {
             a1,
@@ -119,6 +123,7 @@ impl PyNonGravModel {
             m,
             n,
             k,
+            dt,
         })
     }
 
@@ -127,7 +132,7 @@ impl PyNonGravModel {
     ///
     /// See :py:meth:`NonGravModel.new_comet` for more details.
     #[allow(clippy::too_many_arguments)]
-    #[pyo3(signature = (a1, a2, a3, alpha=1.0, r_0=1.0, m= 2.0, n=1.0, k=0.0))]
+    #[pyo3(signature = (a1, a2, a3, alpha=1.0, r_0=1.0, m= 2.0, n=1.0, k=0.0, dt=0.0))]
     #[staticmethod]
     pub fn new_asteroid(
         a1: f64,
@@ -138,6 +143,7 @@ impl PyNonGravModel {
         m: f64,
         n: f64,
         k: f64,
+        dt: f64,
     ) -> Self {
         Self(NonGravModel::JplComet {
             a1,
@@ -148,6 +154,7 @@ impl PyNonGravModel {
             m,
             n,
             k,
+            dt,
         })
     }
 
@@ -166,9 +173,10 @@ impl PyNonGravModel {
                 m,
                 n,
                 k,
+                dt,
             } => format!(
-                "kete.propagation.NonGravModel.new_comet(a1={:?}, a2={:?}, a3={:?}, alpha={:?}, r_0={:?}, m={:?}, n={:?}, k={:?})",
-                a1, a2, a3, alpha, r_0, m, n, k,
+                "kete.propagation.NonGravModel.new_comet(a1={:?}, a2={:?}, a3={:?}, alpha={:?}, r_0={:?}, m={:?}, n={:?}, k={:?}, dt={:?})",
+                a1, a2, a3, alpha, r_0, m, n, k, dt
             ),
         }
     }

src/kete_core/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "kete_core"
-version = "1.0.1"
+version = "1.0.2"
 edition = "2021"
 readme = "README.md"
 license = "BSD-3-Clause"

src/kete_core/src/propagation/nongrav.rs

@@ -1,7 +1,10 @@
 use nalgebra::Vector3;
+use pathfinding::num_traits::Zero;
 
 use crate::constants::{C_AU_PER_DAY_INV_SQUARED, GMS};
 
+use super::analytic_2_body;
+
 /// Non-Gravitational models.
 /// These are used during integration to model non-gravitational forces on particles in
 /// the solar system.
@@ -39,6 +42,9 @@ pub enum NonGravModel {
         n: f64,
         /// Coefficients for the g(r) function defined above.
         k: f64,
+        /// Time delay for the forces, this is applied by propagating the object to the
+        /// specified delay before computing the forces.
+        dt: f64,
     },
 
     /// Dust model, including Solar Radiation Pressure (SRP) and the Poynting-Robertson
@@ -67,6 +73,7 @@ impl NonGravModel {
         m: f64,
         n: f64,
         k: f64,
+        dt: f64,
     ) -> Self {
         Self::JplComet {
             a1,
@@ -77,6 +84,7 @@ impl NonGravModel {
             m,
             n,
             k,
+            dt,
         }
     }
 
@@ -96,6 +104,7 @@ impl NonGravModel {
             m: 2.15,
             n: 5.093,
             k: 4.6142,
+            dt: 0.0,
         }
     }
 
@@ -128,12 +137,18 @@ impl NonGravModel {
                 m,
                 n,
                 k,
+                dt,
             } => {
+                let mut pos = *pos;
                 let pos_norm = pos.normalize();
-                let rr0 = pos.norm() / r_0;
-                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();
+
+                if !dt.is_zero() {
+                    (pos, _) = analytic_2_body(-dt, &pos, vel, None).unwrap();
+                };
+                let rr0 = pos.norm() / r_0;
+                let scale = alpha * rr0.powf(-m) * (1.0 + rr0.powf(*n)).powf(-k);
                 *accel += pos_norm * (scale * a1);
                 *accel += t_vec * (scale * a2);
                 *accel += n_vec * (scale * a3);