Propagation loop never exits for some orbits

[gregjesl]

Bug report

Describe the bug

Propagation loop never exits for some orbits.

To Reproduce

let frame = almanac.frame_from_uid(MOON_PA_FRAME).expect("Could not load frame");
let mut jggrx_meta = MetaFile {
        uri: "https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/jggrx_0660b_sha.tab".to_string(),
        crc32: None
    };
jggrx_meta.process(true)?;
let sph_harmonics = Harmonics::from_stor(
        almanac.frame_from_uid(frame)?,
        HarmonicsMem::from_shadr(&jggrx_meta.uri, 10, 10, false)?,
    );

let epoch = Epoch::from_gregorian_utc_hms(2024, 1, 1, 0, 0, 0);
let orbit = Orbit::try_keplerian_altitude(
        50.0, 
        1e-6, 
        40.5, 
        45.0, 
        0.0, 
        0.0, 
        epoch, 
        frame).unwrap();

let sc = Spacecraft::builder()
            .orbit(orbit)
            .build();

let mut orbital_dyn = OrbitalDynamics::point_masses(vec![]);
orbital_dyn.accel_models.push(harmonics.clone());
let dynamics = SpacecraftDynamics::from_models(orbital_dyn, vec![]);

let future = epoch + Unit::Day * 90;

let (_, trajectory) = Propagator::default(dynamics)
        .with(sc, almanac.clone())
        .until_epoch_with_traj(future)
        .unwrap();

Expected behavior

until_epoch_with_traj should exit after a few seconds.

Platform

Windows 11, PowerShell

Additional context

Another orbit I encountered the issue with was the same parameters but an inclination of 13.500000000000002 and a RAAN of 9.0. Interestingly enough, I did not encounter the issue with an inclination of 13.5.

The orbit is expected to become increasingly eccentric. My initial suspicion is that propagating orbits using harmonics when the orbit falls below the mean radius of the object produces unexpected results.

[gregjesl]

Changing

let future = epoch + Unit::Day * 90;

to

let future = epoch + Unit::Day * 30;

results in the propagator exiting properly. After 30 days the spacecraft would have impacted the moon, so it looks like the issue is only with unrealistic orbits.

In my simulation, I try to catch this case with

until_event(future - epoch, &Event::new(nyx::md::StateParameter::Rmag, 1737.4))

but until_event does exit when the event occurs, it exists once the entire duration is computed. I don't see a method in PropInstance that exits once a condition is reached. I can work around this by implementing a loop that takes smaller timesteps (maybe one day) and checks if the orbit is still valid.

[ChristopherRabotin]

Hi Gregory,

Thanks for the bug report. Have you tried initializing the orbit in the Moon J2000 frame, but keeping the gravity field in the Moon PA frame? The frame of the orbit determines the integration frame, and I don't think the Moon PA frame is an inertial frame so I'm not sure how it would behave.

Concerning the until_event approach, what's the radial magnitude when it exits? And does it take a while to complete? Normally this function should exit properly only if the event is found, otherwise it'll return an error.

I'll recreate your example locally.

[ChristopherRabotin]

This is definitely a bug because the step size is allowed to reach ... zero nanoseconds. I just added the step size in the info log.

 INFO  nyx_space::propagators::instance   > Propagating for 90 days until 2024-03-31T00:00:00 UTC
 INFO  nyx_space::propagators::instance   >     ... current epoch 2024-03-16T03:39:12.541930440 UTC, remaining 14 days 20 h 20 min 47 s (step size = 0 ns)

I'll fix this right away, thanks for the report.

[ChristopherRabotin]

In the following code, I built upon your example and used the J2000 frame for propagation, and kept the harmonics in the Moon PA frame (as they should be). I also used your approach to propagate until impact on the Moon using the event finder. It works fine, but I'll still push the step size fix in a moment.

_Note: I've built this on top of the LRO example in my local clone of the repo, so there are plenty of unused imports. I'm also using the Moon BPC moon_pa_de440_200625.bpc which doesn't have MOON_PAFRAME as ID 31000 but instead as 31008, so I've done the quick fix for that. Without the fix, Nyx errors saying it doesn't know how to rotate from Moon J2000 (frame in which the spacecraft is defined) into the Moon Principal Axes frame, and that transformation is required to compute the spherical harmonics in the other frame

╰─⠠⠵ RUST_LOG=info cargo run --example gh377 --release
(...)
     Running `target/release/examples/gh377`
 INFO  anise::almanac::metaload::metafile > Using cached /home/chris/.local/share/nyx-space/anise/de421.bsp
 INFO  anise::almanac::metaload::metafile > Using cached /home/chris/.local/share/nyx-space/anise/pck11.pca
 INFO  anise::almanac::metaload::metafile > Using cached /home/chris/.local/share/nyx-space/anise/moon_pa_de440_200625.bpc
 INFO  anise::almanac::metaload::metafile > Using cached /home/chris/.local/share/nyx-space/anise/lrorg_2023349_2024075_v01_LE.bsp
 INFO  anise::almanac                     > Loading almanac from /home/chris/.local/share/nyx-space/anise/de421.bsp
 INFO  anise::almanac                     > Loading as DAF/SPK
 INFO  anise::almanac                     > Loading almanac from /home/chris/.local/share/nyx-space/anise/pck11.pca
 INFO  anise::almanac                     > Loading almanac from /home/chris/.local/share/nyx-space/anise/moon_pa_de440_200625.bpc
 INFO  anise::almanac                     > Loading as DAF/PCK
 INFO  anise::almanac                     > Loading almanac from /home/chris/.local/share/nyx-space/anise/lrorg_2023349_2024075_v01_LE.bsp
 INFO  anise::almanac                     > Loading as DAF/SPK
 INFO  anise::almanac::metaload::metafile > Using cached /home/chris/.local/share/nyx-space/anise/jggrx_0660b_sha.tab
 INFO  nyx_space::io::gravity             > /home/chris/.local/share/nyx-space/anise/jggrx_0660b_sha.tab loaded with (degree, order) = (10, 10)
 INFO  nyx_space::propagators::instance   > Searching for Rmag = 1.7374e3 km (± 1e-3 km)
 INFO  nyx_space::propagators::instance   > Propagating for 90 days until 2024-03-31T00:00:00 UTC
 INFO  nyx_space::propagators::instance   >     ... done in 30 s 862 ms 492 μs 115 ns
 INFO  nyx_space::md::events::search      > Searching for Rmag = 1.7374e3 km (± 1e-3 km) with initial heuristic of 21 h 36 min
 INFO  nyx_space::md::events::search      > Event Rmag = 1.7374e3 km (± 1e-3 km) found 65 times from 2024-01-16T04:54:13.338489790 UTC until 2024-03-30T02:57:26.422865879 UTC
total mass = 0.000 kg @  [Moon J2000] 2024-01-16T04:54:13.338489790 UTC sma = 1787.676999 km    ecc = 0.029394  inc = 45.163035 deg     raan = 34.241078 deg    aop = 263.634316 deg    ta = 17.403958 deg  Coast
RMAG = 1737.400505842332 km

Here's the code that worked:

extern crate nyx_space as nyx;
extern crate pretty_env_logger as pel;

use anise::{
    almanac::metaload::MetaFile,
    constants::{
        celestial_objects::{EARTH, JUPITER_BARYCENTER, MOON, SUN},
        frames::{EARTH_J2000, MOON_J2000, MOON_PA_FRAME},
    },
};
use hifitime::{Epoch, TimeUnits, Unit};
use nyx::{
    cosmic::{Aberration, Frame, MetaAlmanac, SrpConfig},
    dynamics::{
        guidance::LocalFrame, Harmonics, OrbitalDynamics, SolarPressure, SpacecraftDynamics,
    },
    io::{ConfigRepr, ExportCfg},
    md::{
        prelude::{HarmonicsMem, Traj},
        Event,
    },
    od::{
        msr::RangeDoppler,
        prelude::{TrackingArcSim, TrkConfig, KF},
        process::{Estimate, NavSolution, ODProcess, ResidRejectCrit, SpacecraftUncertainty},
        snc::SNC3,
        GroundStation,
    },
    propagators::Propagator,
    Orbit, Spacecraft, State,
};

use std::{collections::BTreeMap, error::Error, path::PathBuf, str::FromStr, sync::Arc};

fn main() -> Result<(), Box<dyn Error>> {
    pel::init();

    let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
        .iter()
        .collect();

    let meta = data_folder.join("lro-dynamics.dhall");

    // Load this ephem in the general Almanac we're using for this analysis.
    let mut almanac = Arc::new(
        MetaAlmanac::new(meta.to_string_lossy().to_string())
            .map_err(Box::new)?
            .process(true)
            .map_err(Box::new)?,
    );

    let frame = almanac
        .frame_from_uid(MOON_PA_FRAME.with_orient(31008))
        .expect("Could not load frame");
    let mut jggrx_meta = MetaFile {
        uri: "https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/shadr/jggrx_0660b_sha.tab".to_string(),
        crc32: Some(0xebdb7283)
    };
    jggrx_meta.process(true)?;
    let sph_harmonics = Harmonics::from_stor(
        almanac.frame_from_uid(frame)?,
        HarmonicsMem::from_shadr(&jggrx_meta.uri, 10, 10, false)?,
    );

    let epoch = Epoch::from_gregorian_utc_hms(2024, 1, 1, 0, 0, 0);
    let orbit = Orbit::try_keplerian_altitude(
        50.0,
        1e-6,
        40.5,
        45.0,
        0.0,
        0.0,
        epoch,
        almanac.frame_from_uid(MOON_J2000).unwrap(),
    )
    .unwrap();

    let sc = Spacecraft::builder().orbit(orbit).build();

    let mut orbital_dyn = OrbitalDynamics::point_masses(vec![]);
    orbital_dyn.accel_models.push(sph_harmonics.clone());
    let dynamics = SpacecraftDynamics::from_models(orbital_dyn, vec![]);

    let future = epoch + Unit::Day * 90;

    let (sc, trajectory) = Propagator::default(dynamics)
        .with(sc, almanac.clone())
        .until_event(
            future - epoch,
            &Event::new(nyx::md::StateParameter::Rmag, 1737.4),
        )
        .unwrap();
    // .until_epoch_with_traj(future)
    // .unwrap();

    println!("{sc:x}");
    println!("RMAG = {} km", sc.orbit.rmag_km());

    Ok(())
}

I'll open a PR momentarily with the fix.