Forming a Nyx::Orbit from Keplerian parameters

[gwbres]

Hello @ChristopherRabotin,

RINEX Navigation files (so called BRDC NAV Files for broadcast navigation ) are a snapshot of radio messages streamed by satellite vehicles at some point in time. They are used to later on determine the SV position, and we need 4 of them to resolve a position on the ground. The most fundamental frame contains the Keplerian parameters of the SV orbit at this point in time, and we apply the Kepler equation to determine {x, y, z}.

To do so, I implemented the Kepler solver that ESA recommends, we use Ephemeris::kepler2ecef(t) that returns {x, y, z} at t .

I would like to get rid of kepler2ecef and use your API to do the work

• it would reduce the code quantity
• it has been a difficult function to develop & very difficult to verify
• I suspect it does not work fine for "exotic" constellations. Relying on your library (also thorougly verified), will most likely fix that
• I'm already able to form a nyx::Orbit in other parts of my ecosystem where {x, y, z} have already been determined, and then take advantage Nyx for more advanced stuff (like Eclipse conditions determination)

But the problem is I don't understand the way you formulate keplerian parameters, as it does not look like the parameters we typicall get over the radio. Could you help me figure out the translation ?

Keplerian parameters

These are the typical parameters the satellites stream over the radio. Many of them, I am unsure about their meaning

• sqrt(a): square root of the semi major axis [sqrt(m)]
• Delta n [radians / sec]: pretty sure this is d(mean motion)/dt
• M0 [radians]: not documented. I think this is the mean anomaly but I'm far from sure.
  Consider that the solver implements (which you may recognize):

        // `a` is sqrt(semi major axis) 
        // `dt`  is `t` the Epoch we would like to resolve {x, y, z} expressed as seconds within its timescale
        let n0 = (EARTH_GRAVITATIONNAL_CONSTANT / kepler.a^3).sqrt()
        let n = n0 + delta_n;
        let omega_k = kepler.m0 + n *  dt; // << M0

• capital omega0 [radians] and capital omega dot in [radians/sec] (dot=derivative in time).
  Not documented. It is not clear whether this is the longitude of the ascending node, or the argument of periapsis.
  The capital makes me lean towards the first.
  Consider that the solver implements (which you may recognize):

        // `TOE` is the reference point in Time (epoch) for the current Ephemeris period
        // `dt`  is `t` the Epoch we would like to resolve {x, y, z} expressed as seconds within its timescale
           let omega_k = kepler.omega_0
            + (perturbations.omega_dot - Kepler::EARTH_OMEGA_E_WGS84) * dt
            - Kepler::EARTH_OMEGA_E_WGS84 * kepler.toe;

• i0 and d(i0)/dt [radians] these are the only terms I'm certain what they are: inclination and rate of inclination angles.
• Osculating terms, expressed as one term in meters, one cosine and one sine term both in radians,
  that we implement on all three axis
• tiny omega, which I think is the argument of periapsis, and is badly named. Expressed in radians, that the recommended calculations uses right there

That's it. If you teach me how to convert those into the following parameters, we can get rid of one complex task

see you 🚀

[gwbres]

from the specs, I'm figuring:

• dt: obvious
• Frame: you taught me last year
• sma: semi major axis, therefore a^2 in my case
• aop: argument of periapsis, which I think is tiny omega in RINEX notations
• ta: true anomaly which I need to compute from mean anomaly (M0) and I have not figured how to do that yet
• raan: seems like capital omega in RINEX notation
• inc: not sure whether it is i0 directly or i_0 + i_dot * dt + osculating/perturbation(di)

Guess I'll run a couple of tests..

[ChristopherRabotin]

Hi Guillaume,

If you can, I would recommend using ANISE because the Earth rotation model is sure to be correct. In fact, if you have historical data, you could even use the ITRF93 Earth body fixed frame, which is attached to the crust of the Earth and computed from the latest observations (published daily by NASA). If you need the features of Nyx for your computations (like the filter or the propagation), and not just what's in ANISE (frame transformations, Keplerian computations, etc.), then I would recommend that you (temporarily) use the https://github.com/nyx-space/nyx/tree/86-use-spice-files-via-anise-for-importexport-and-computations branch directly. I hope to finish that chunk of work soon.

Back to your question. Whether you're using the latest beta of Nyx without ANISE, or ANISE itself, you'll want to use the try_keplerian_mean_anomaly initialization of an Orbit. From my understanding of page 21 of this pdf, the (capital) Omega is the right ascension of the ascending node (i.e. the angle difference between the "zero angle" (vernal equinox pointing) and where the spacecraft crosses the equatorial plane). The small omega is the argument of periapse (which can be seen as the "orientation" of the orbit's axis).

So here is what I recommend (I'm making the assumption that you're using either "pure" ANISE or the Nyx 86-blah-blah branch):

1. Load the latest Almanac using MetaAlmanac::latest() . This will fetch the very latest Earth orientation parameters, along with the latest planetary constants (mass and shape of Earth, etc.).
2. Grab the latest high precision Earth body fixed frame, ITRF93 (or you can use the lower precision frame called "IAU Earth") using frame_from_uid: let earth_hp = almanac.frame_from_uid(EARTH_ITRF93).expect("frame not available");
3. Compute the mean anomaly at the specified time using this part of your code: https://github.com/georust/rinex/blob/main/rinex/src/navigation/ephemeris.rs#L364 . Compute the RAAN at the specified time using the Omega_0 and the Earth Omega drift rate in rad/s (be sure to convert all angles from radians to degrees for the Orbit initialization and the SMA from meters to kilometers).
4. And finally, build the Orbit structure: let orbit = Orbit::try_keplerian_mean_anomaly(sma_km, ecc, inc_deg, raan_deg, aop_deg, ma_deg, epoch, earth_hp).expect("initialization failed?!"); // earth_hp is the Earth high precision frame

Additional notes:

• Since everything you're using is on hifitime v4, you would need to use ANISE's latest master branch from this morning.
• If you're using the 86-... branch of Nyx, note that the paths in the Cargo.toml are relative and don't use the git paths! This is because it helps my debugging of that branch, but all of the paths will be fixed before the merge. But if you're using that branch, you'll need to update the Cargo.toml branch too.

Let me know how this goes.

[gwbres]

If you could do that, that would be great indeed.
I was just checking, I'm still not finished working on other parts of my crates, I won't have time to investigate in the next 10 days I would say

[ChristopherRabotin]

I'll work on it this weekend, I'm a bit overloaded until then.

[gwbres]

Hello Chris, no worries, I'm still working on the improvements of our RINEX/SP3 analysis mode, before getting back to positioning

[ChristopherRabotin]

Oof, finally! georust/rinex#252 Sorry for the delay!

[gwbres]

I was not expecting you'd contribute on this topic 🚀
thank's a lot, it will surely help the process, follow up in the new discussion.
I will close this one when we select a solution