GMAT vs nyx difference in output states

[skincell]

I'm a novice in Rust, so please bear with me. I ran into a difference between a propagated state using only Earth point mass model. The nyx propagator that I used was CashKarp45 and the RK89 from GMAT. This could be due to the difference in the propagators rather than a nyx error; however, I put this here for documentation purposes just in case. The options that are specified for the propagator initialization is the following:

Propagator::new::<CashKarp45>(&Options::with_adaptive_step(0.1, 30.0, 1e-12));

Input PropInfo { init: Some(State { datetime: 58112.04877314817, x: 6678.0, y: 0.0, z: 0.0, x_dot: -1.034425773665938, y_dot: 8.470618983406897, z_dot: 0.0 }), gm: 398600.4415, prop_time: 300.0 } Output Sending: State { datetime: 58112.05224537039, x: 5959.370605933568, y: 2488.2135004453403, z: 0.0, x_dot: -3.7494156073202514, y_dot: 7.926583218577195, z_dot: 0.0 }

The script is attached. You will have to rename it as a Test.script file. I put it as a .txt file to upload it here. Test.txt

The output state for the 300 second GMAT propagation was: X: 5959.745504208833, Y: 2487.420821985822, Z:0 VX:-3.748533621228721, VY: 7.926951403625875, VZ:0 Elapsed Time: 299.9999997206032

The Y coordinate differs by almost 1 km for a 300 second propagation. The good news is that when I switched my propagator to a RK87 using nyx. The results differed at the sub-meter level.

I am using nyx-space 0.0.3, so this issue might be resolved in the latest builds.

[ChristopherRabotin]

It seems like you are using a Prince Dormand 45 in GMAT but a Cash Karp 45 in nyx. These are different integrators, and the results will absolutely vary between them. Can you check that by using the Dormand45 integrator in nyx you get the same output state that GMAT returns? To do this, change the "use" instruction which imports the CashKarp45 to add Dormand45, and change the line in the code where CashKarp45 is used.

Christopher Rabotin.

On Tue, Apr 23, 2019, 13:14 skincell notifications@github.com wrote:

I'm a novice in Rust, so please bear with me. I ran into a difference between a propagated state using only Earth point mass model. The nyx propagator that I used was CashKarp45 and the RK89 from GMAT. This could be due to the difference in the propagators rather than a nyx error; however, I put this here for documentation purposes just in case. The options that are specified for the propagator initialization is the following:

Propagator::new::(&Options::with_adaptive_step(0.1, 30.0, 1e-12));

Input PropInfo { init: Some(State { datetime: 58112.04877314817, x: 6678.0, y: 0.0, z: 0.0, x_dot: -1.034425773665938, y_dot: 8.470618983406897, z_dot: 0.0 }), gm: 398600.4415, prop_time: 300.0 } Output Sending: State { datetime: 58112.05224537039, x: 5959.370605933568, y: 2488.2135004453403, z: 0.0, x_dot: -3.7494156073202514, y_dot: 7.926583218577195, z_dot: 0.0 }

The script is attached. You will have to rename it as a Test.script file. I put it as a .txt file to upload it here. Test.txt https://github.com/ChristopherRabotin/nyx/files/3108975/Test.txt

The output state for the 300 second GMAT propagation was: X: 5959.745504208833, Y: 2487.420821985822, Z:0 VX:-3.748533621228721, VY: 7.926951403625875, VZ:0 Elapsed Time: 299.9999997206032

The Y coordinate differs by almost 1 km for a 300 second propagation. The good news is that when I switched my propagator to a RK87 using nyx. The results differed at the sub-meter level.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/ChristopherRabotin/nyx/issues/58, or mute the thread https://github.com/notifications/unsubscribe-auth/ABEZV2BEZM5MBOVPDG3SXK3PR5GXPANCNFSM4HH5E5RA .

[ChristopherRabotin]

@skincell, any update on the test I mention above?

[skincell]

I have updated to nyx 0.5, and here is the output for the same input state as before using the Dormand45: adaptive step: 0.001 minimum step, 60 largest step, and 1e-15 tolerance. end state after propagating for 300 seconds.

Input DEBUG goats > Received: PropInfo { init: Some(State { datetime: 58112.04877314817, x: 6678.0, y: 0.0, z: 0.0, x_dot: -1.0343754813787296, y_dot: 8.470617138524924, z_dot: 0.0 }), gm: 398603.1, prop_time: 300.0 } Output DEBUG goats > Sending: State { datetime: 58112.05224537039, x: 5958.859115788939, y: 2489.3211942845674, z: 0.0, x_dot: -3.7506116399140694, y_dot: 7.9260649205746025, z_dot: 0.0 }

Output for GMAT as previous: X: 5959.745504208833, Y: 2487.420821985822, Z:0 VX:-3.748533621228721, VY: 7.926951403625875, VZ:0 Elapsed Time: 299.9999997206032

---

I also ran into issues with the RK89 not matching the output from GMAT's RK89 at the 100s of meters level. I ended up changing it to RK4 fixed for nyx to match the RK89 better than the meter level which is what I'm looking for. The input state that I see this difference for is:

DEBUG goats > Received: PropInfo { init: Some(State { datetime: 58112.04877314817, x: 6677.97, y: 0.0, z: 0.0, x_dot: -1.0343579648824441, y_dot: 8.47064538331271, z_dot: 0.0 }), gm: 398603.1, prop_time: 1200}.

Let me know if you want the inputs and outputs that I'm getting on my side for this second case (if you think it is pertinent to the issue.)

[ChristopherRabotin]

I'll have to write up a quick demo of this bug. I'm quite surprised because there are validation test cases between nyx and GMAT which show that the accuracy is down to the micrometer.

Christopher Rabotin.

On Thu, May 2, 2019, 10:01 skincell notifications@github.com wrote:

I have updated to nyx 0.5, and here is the output for the same input state as before using the Dormand45: adaptive step: 0.001 minimum step, 60 largest step, and 1e-15 tolerance. end state after propagating for 300 seconds.

Input DEBUG goats > Received: PropInfo { init: Some(State { datetime: 58112.04877314817, x: 6678.0, y: 0.0, z: 0.0, x_dot: -1.0343754813787296, y_dot: 8.470617138524924, z_dot: 0.0 }), gm: 398603.1, prop_time: 300.0 } Output DEBUG goats > Sending: State { datetime: 58112.05224537039, x: 5958.859115788939, y: 2489.3211942845674, z: 0.0, x_dot: -3.7506116399140694, y_dot: 7.9260649205746025, z_dot: 0.0 }

Output for GMAT as previous: X: 5959.745504208833, Y: 2487.420821985822, Z:0 VX:-3.748533621228721, VY: 7.926951403625875, VZ:0 Elapsed Time: 299.9999997206032

I also ran into issues with the RK89 not matching the output from GMAT's RK89 at the 100s of meters level. I ended up changing it to RK4 fixed for nyx to match the RK89 better than the meter level which is what I'm looking for. The input state that I see this difference for is:

DEBUG goats > Received: PropInfo { init: Some(State { datetime: 58112.04877314817, x: 6677.97, y: 0.0, z: 0.0, x_dot: -1.0343579648824441, y_dot: 8.47064538331271, z_dot: 0.0 }), gm: 398603.1, prop_time: 1200}.

Let me know if you want the inputs and outputs that I'm getting on my side for this second case (if you think it is pertinent to the issue.)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/ChristopherRabotin/nyx/issues/58#issuecomment-488731444, or mute the thread https://github.com/notifications/unsubscribe-auth/ABEZV2CDHCN2NIAAV3VSKNTPTMF7LANCNFSM4HH5E5RA .

[skincell]

Let me know if you can reproduce it. It could totally be on my side as I'm unfamiliar with rust. Thanks!

[ChristopherRabotin]

DEBUG goats > Received: PropInfo { init: Some(State { datetime: 58112.04877314817, x: 6677.97, y: 0.0, z: 0.0, x_dot: -1.0343579648824441, y_dot: 8.47064538331271, z_dot: 0.0 }), gm: 398603.1, prop_time: 1200}. ==> notice how the GM you are using is not the default one of GMAT, nor the default one from nyx. The default is:

[src/main.rs:27] EARTH::gm() = 398600.4415

I've coded up a script which shows that GMAT and nyx are within 1e-6 km of difference in radius and 1e-9 km/s in velocity when using a DormandPrince 45 propagator.

Output

cargo run
   Compiling nyxval v0.1.0 (/home/chris/Workspace/DSM/nyxval)
    Finished dev [unoptimized + debuginfo] target(s) in 1.79s
     Running `target/debug/nyxval`
[src/main.rs:27] EARTH::gm() = 398600.4415
Radius error: 5.398898e-7 km
Velocity error: 5.889294e-10 km/s
IntegrationDetails { step: 0.4161779199746647, error: 0.000000000000595205456465188, attempts: 1 }

Reproducibility:

1. Create a folder, e.g. mkdir nyxval
2. Create a file called Cargo.toml and paste the content from the section below
3. Create a folder called src: mkdir src
4. Copy the contents of the main.rs file in there
5. cargo run

Cargo.toml

[package]
name = "nyxval"
version = "0.1.0"
authors = ["Christopher Rabotin <rabotin@advanced-space.com>"]
edition = "2018"

[dependencies]
nalgebra = "0.17"
nyx-space = { git = "https://github.com/ChristopherRabotin/nyx", branch="master"}

main.rs

extern crate nalgebra as na;
extern crate nyx_space as nyx;

use na::{Vector6, U3};
use nyx::celestia::{CelestialBody, EARTH};
use nyx::dynamics::celestial::TwoBody;
use nyx::propagators::error_ctrl::RSSStatePV;
use nyx::propagators::*;

fn main() {
    let init =
        Vector6::from_row_slice(&[6678.0, 0.0, 0.0, -1.034425773665938, 8.470618983406897, 0.0]);
    let rslt = Vector6::from_row_slice(&[
        5959.745504208833,
        2487.420821985822,
        0.0,
        -3.748533621228721,
        7.926951403625875,
        0.0,
    ]);

    let prop_time = 299.9999997206032;
    let min_step = 0.001;
    let max_step = 2700.0;
    let accuracy = 1e-12;

    dbg!(EARTH::gm());

    let mut dynamics = TwoBody::from_state_vec::<EARTH>(init.clone());
    let mut prop = Propagator::new::<Dormand45>(
        &mut dynamics,
        &PropOpts::with_adaptive_step(min_step, max_step, accuracy, RSSStatePV {}),
    );
    prop.until_time_elapsed(prop_time);
    let final_state = prop.state().clone();

    let err_radius = (&final_state.fixed_rows::<U3>(0).into_owned()
        - &rslt.fixed_rows::<U3>(0).into_owned())
        .norm();

    let err_velocity = (&final_state.fixed_rows::<U3>(3).into_owned()
        - &rslt.fixed_rows::<U3>(3).into_owned())
        .norm();

    println!(
        "Radius error: {:.6e} km\nVelocity error: {:.6e} km/s",
        err_radius, err_velocity
    );

    println!("{:?}", prop.latest_details());
}