Fix the RINEX/GNSS tools

[gwbres]

Hello,

the RINEX/GNSS tools are broken since upgrade to the new ANISE+NYX combination (about one month), we need some serious help.
The impossibility to combine hifitime-4 to previous NYX only has caused our previous recent tags to break.
We are facing three problems. Let's fix the main branch and I will apply the corrections to other topics

Precise navigation

The final problem we need to solve is that each epoch (when navigating) terminates in the new Error::Physics nav-solutions/gnss-rtk@80414f3.

In all navigation scenarios.
But it is simply the consequence of the following "problems".

RINEX only navigation with ANISE

The RINEX navigation core fails to resolve positions. When navigating without SP3, we then end up in Error::Physics.

nav-solutions/rinex@7fd3ac5

git clone https://github.com/georust/rinex
cargo build --release --all-features

./target/release/rinex-cli \
-f test_resources/CRNX/V3/MOJN00DNK_R_20201770000_01D_30S_MO.crx.gz\
 -f test_resources/NAV/V3/MOJN00DNK_R_20201770000_01D_MN.rnx.gz \
-f test_resources/CLK/V3/GRG0MGXFIN_20201770000_01D_30S_CLK.CLK.gz \
-P GAL \
-P C1C,C5Q,L1C,L5Q -p

[2024-07-20T09:23:06Z DEBUG gnss_rtk::solver] 2020-06-25T07:02:30 GPST (E36) : InterpolationResult { elevation: 0.0, azimuth: 0.0, position: [[27116923.03778262, -11859500.49457391, 10167.250998856873]], velocity: None }
[2024-07-20T09:23:06Z WARN  rinex_cli::positioning::ppp] 2020-06-25T07:02:30 GPST : pvt solver error "physics issue: retrieving gravitational parameter requires the frame Earth J2000 to have mu_km3_s2 defined"

To be more accurate, the ECEF position of the SV has been output, it is most likely incorrect (unit/scaling error ?), this is keplerian navigation. The fact we're stuck at elevation=azmith=0 means that the proposal is incorrect and the solver never achieves initialization.
We remain in that state for ever and never generate a single output.

SP3 navigation with ANISE

Precise navigation is not able to deploy because I'm unable to use the latest API, or more accurately, have no idea what model we need to load and how we should do that

git clone https://github.com/georust/rinex
cargo build --release --all-features

./target/release/rinex-cli \
-f test_resources/CRNX/V3/MOJN00DNK_R_20201770000_01D_30S_MO.crx.gz\
-f test_resources/NAV/V3/MOJN00DNK_R_20201770000_01D_MN.rnx.gz \
-f test_resources/CLK/V3/GRG0MGXFIN_20201770000_01D_30S_CLK.CLK.gz \
-f test_resources/SP3/V3/GRG0MGXFIN_20201770000_01D_15M_ORB.SP3.gz \
-P GAL \
-P C1C,C5Q,L1C,L5Q -p

[2024-07-20T09:21:19Z TRACE reqwest::blocking::client] closed runtime thread (ThreadId(2))
thread 'main' panicked at rinex-cli/src/positioning/orbit/sp3.rs:79:33:
failed to build Almanac: processing file #0 (MetaFile { uri: "http://public-data.nyxspace.com/anise/v0.4/de440s.bsp", crc32: Some(3072159656) }) caused an error: fetching http://public-data.nyxspace.com/anise/v0.4/de440s.bsp returned 404 Not Found
stack backtrace:

To fix that, the model we download must be corrected and obtained

BeiDou case

BeiDou navigation is very particular. Some vehicles are geostationnary vehicles (it is determined from the Vehicle #ID number).
I suspect the rinex::navigation keplerian modifications to ANISE 4 has broken the support of BeiDou:
nav-solutions/rinex@7fd3ac5

My expectations is that the ANISE 4 upgrade should apply to LEO Beidou and simplify the previous process, but we will have to re-introduce the GEO logic that has been lost in that contribution

[ChristopherRabotin]

Hi Guillaume,

I'm sorry you're having issues since the upgrade. I will help fix them. Which of the repos should I work on first? I hope we can fix this all this weekend.

Best,
Chris

[ChristopherRabotin]

This sounds promising. Where is this ".orbit()" call? In gnss_rtk ? This could be the frame issue!

[gwbres]

I just improved my previous statement

[ChristopherRabotin]

This is very helpful. I can dive into this today.

[gwbres]

if orbit.ecc() is an angle, it's probably missing a conversion to radians ?
The error is in the SPEED_OF_LIGHT_KM_S which is wrong and should be SPEED_OF_LIGHT_M_S
This is also the explnation why dt is so boggous, it should be M_S and not KM_S.

All gnss_rtk world is expressed in meters (SI units) to simplify things (and Fixed body frames).
The principle of the (ground based) navigation equation, is to consider a fixed reference system (=body fixed frame?) on said planet, and solving in a "differential" manner, that means, all results are expressed with respect to an apriori position on the ground.

[ChristopherRabotin]

Ecc is the eccentricity, which does not have any units : https://nyxspace.com/anise/reference/mathspec/orbital_elements/#eccentricity-evec-ecc

[gwbres]

Thank you for helping me out, we're almost done !

The last topic to solve is replacing our rinex:navigation:ephemeris keplerian calculations by the Keplerian solver integrated within ANISE.

Your data analysis on what we forward to Orbit::try_keplerian was correct: the modified quantities and scaling are mandatory (I kept your comments in the code).

The performance is degraded by about 10x/100x times I would say, currently, so radio based navigation is temporarily broken.
Our previous steps have fixed other navigation scenarios.
The error is purely located rinex:navigation:ephemeris, because it is solely due to the SV state vector and other navigation scenarios do work. You can try it out with :

# on rinex (main)
./target/release/rinex-cli \ 
   -f test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz \
   -f test_resources/NAV/V3/ESBC00DNK_R_20201770000_01D_MN.rnx.gz \
   -P GPS \ 
   -p

Our reference point, for radio based navigation, was this contribution :
nav-solutions/rinex#249. We still have the option to backport it, the switch to ANISE in this portion, simplifies what the RINEX lib is reponsible for.

It puzzles me that the resolved orbital states are not very incorrect, in fact, if we run something like this, we can compare the orbital state resolved from Keplerian calculations, to the state parsed in the SP3 file. Sometimes we're way off, but it does not happen this much, and the compass projection looks very fine

./target/release/rinex-cli \ 
   -f test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz \
   -f test_resources/NAV/V3/ESBC00DNK_R_20201770000_01D_MN.rnx.gz \
  -f test_resources/SP3/GRG0MGXFIN_20201770000_01D_15M_ORB.SP3.gz \
   -P GPS \ 
   -g --sky --orbit

[gwbres]

I'm investigating the last remaining error and it is harder than the previous ones.

As a reminder, we used to have performant radio navigation using MEO vehicles, now it's partially broken.
I'm comparing past results to the latest.
If I compare the SV {x, y, z} state vector that is obtained now (from Kepler) and previous:

git clone https://github.com/georust/rinex
cargo build --all-features --release
./target/release/rinex-cli \
   -f test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz \
  -f test_resources/NAV/V3/ESBC00DNK_R_20201770000_01D_MN.rnx.gz \
  -P GPS \
  -p

• we're getting x=17921612, y=5139682, z=18837047 [ECEF_m]
• while we used to obtain x=17921489, y=5139473, z=18837044 [ECEF_m]

We're very close, about 200m error on x and y, and z is almost identical.
Note that we cannot consider the past results as absolute truth either, due to the complexity of the task to obtain correct keplerian navigation from RINEX interpretation. So it's possible that the new z is the truth.

These are the input to gnss_rtk, and any error on them reflects on the following process and degrades the solution.
Note that gnss_rtk expects ECEF coordinates (as of today).

Now I can do this to investigate a little further. We introduce an SP3 and run the Radio/SP3 comparison mode.

./target/release/rinex-cli \
   -f test_resources/CRNX/V3/ESBC00DNK_R_20201770000_01D_30S_MO.crx.gz \
  -f test_resources/NAV/V3/ESBC00DNK_R_20201770000_01D_MN.rnx.gz \
  -P GPS \
  -g --orbit --orbit-residual

The 3D proj does not emphasize anything because we"re too close from truth (previous statement)

The x, y, z comparison between SV { x, y, z } obtained from Kepler expressed in ECEF, and SV { x , y , z } parsed from SP3 in slightly different frames (IGb14) shows outstanding results, which tends to demonstrate that the calculation is correct

So right now I am kind of puzzled.
There is either an error in the Time instant considered to express the state vector, or there is an error introduced by the Almanac knowledge within gnss_rtk and the state vector we provide here

[gwbres]

Do you think it would be possible to integrate the derivative terms, which must be considered according to the kepler calculation formula and the RINEX specs, and that we do not take into account in the new try_keplerian method ? the derivative terms are calculated right new to it

There were also some tiny rotations that we used to apply when converting to ECEF frame specifically:

[gwbres]

actually, if I backport the previous correct state, and run the comparison plot,
i'm getting 100x times better results !

[gwbres]

I have backported to the latest stable state and left the anise port commented out inside this code

[ChristopherRabotin]

I'm not understanding the plots correctly because I see an error on the order of a micrometer, which looks pretty good to me.

The rotation in the meo_orbit_to_ecef looks like it's the Hill frame rotation, which is a 3-1-3 rotation by the axes of the right ascension of the ascending node on the Z axis, then a rotation by the inclination angle on the X axis, and then usually a final rotation by Z corresponding to the phasing angle used (.e.g either true anomaly for non-circular orbits, or true longitude otherwise). This frame is usually used for a quick approximation of the conversion from ECEF to ECI, but it isn't high fidelity. ANISE should provide high fidelity rotations, so I'm quite confused that you're seeing better results with a lower fidelity rotation than with a higher fidelity one.

In this commit, nav-solutions/rinex@0fae4d8, it seems that you've entirely removed the use of ANISE. Is that right? My guess is that there's a frame difference somewhere where one of the functions expects a given frame but receives a position in another frame, and the rotation from the meo_orbit function masks that away.

With ANISE, when using the daily publication of the high precision frame, cf. https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/earth_latest_high_prec.cmt, you should get essentially zero error compared to the SP3 data because the high precision frame is built from the daily SP3 data of dozens of ground stations around the world. This is the objective I would like to achieve when combining RINEX + ANISE.

[gwbres]

I'm not understanding the plots correctly because I see an error on the order of a micrometer, which looks pretty good to me.

It is supposed to be the difference between the mass center of the spacecraft obtained by kepler nav, and the one read from the SP3 file. Now that I've said that, it does not make sense, from both outstanding performances, that one diverges this much.
Although we do have "worse" results from the last ANISE attempt

The rotation in the meo_orbit_to_ecef looks like it's the Hill frame rotation, which is a 3-1-3 rotation by the axes of the right ascension of
the ascending node on the Z axis, then a rotation by the inclination angle on the X axis, and then usually a
final rotation by Z corresponding to the phasing angle used (.e.g either true anomaly for non-circular orbits,
or true longitude otherwise). This frame is usually used for a quick approximation of the conversion from ECEF to ECI,
but it isn't high fidelity. ANISE should provide high fidelity rotations, so I'm quite confused that you're seeing better
results with a lower fidelity rotation than with a higher fidelity one.

Very interesting. At first I only saw ANISE as a simple code simplifier. But it is clear that it should also come with performance enhancement. One key element being the high fidelity rotations.
Although for the Keplerian calculator it does not simplify much, due to the little integration of the raw RINEX data fields we have to do (RAAN and other fields).

In this commit, nav-solutions/rinex@0fae4d8, it seems that you've entirely removed the use of ANISE. Is that right? My guess is that there's a frame difference somewhere where one of the functions expects a given frame but receives a position in another frame

I simply commented it out, so most of your work and data investigation is not lost and can be brought back quickly.
It has taught me how to deploy and interface to the basics of ANISE, which is already important.
My strategy right now is to investigate how to get rid of my Elevation / Azimuth calculator and replace it by ANISE. It could be a first step, and is much easier than the keplerian navigation. We don't have much rotations other than that, and all our calculations are ECEF based. Are we sure we can rule out a bad scaling/unit error yet ?

[gwbres]

The ANISE port is almost done I think, it's giving correct results now.

I created a new branch (orbit) where I first worked on the new orbital projection reporting (check them out).
Then I gave the ANISE import a new attempt, and it seems to be working fine.

You can try 4 things. Unfortunately the reporting logic is not 100 % correct, depending on the Input combination and/or filter designer, some plots are missing. If you run those, you will see all the possible plots at least once

git clone https://github.com/georust/rinex -b orbit
cargo build --all-features -r

# Static/BRDC/Galileo  x_err=-128mm, y_err=71mm, z_err=40mm
# Nothing new here, outstanding results obtained 2 months ago when BRDC nav was finally fixed.
# Now it's using ANISE and the performance is preserved
./tutorials/GAL/esberg-brdc.sh

# Static/BRDC/GPS  x_err=-56cm, y_err=70cm, z_err=-1.2m
# Best results obtained ever with GPS (until today we could not get below 1m)
./tutorials/GPS/esbjerg-brdc.sh

# This SP3 analysis will unlock the orbit projections.
# In this case, the orbit is parsed from the text file, attitudes are resolved with Anise, elevation/angle comes from
# rinex:navigation:ephemeris, which uses ANISE for that purpose now. 
# The results are wrong since we're able to obtain negative elevation.
# If that is the true behavior, that it will degrade the navigation performances, since elevation is used in physical/environmental perturbations modeling.
./tutorials/SP3/basic.sh

# This SP3+BRDC analysis will unlock the residual error plot - that I corrected today.
# The error is calculated by reading the high precision coordinates from SP3, solving Kepler from radio Ephemeris 
# at that very epoch, for the same vehicle, and substract the two.
./tutorials/SP3/basic+brdc.sh

I also uses the Orbit.lat/lon projection directly from ANISE now, and not from map_3d.

Since BRDC navigation has been fixed (a few weeks ago), I'm very surprized by how well it performs compared to SP3.
It surpassed by a factor of x10. While it should actually be the other way around.
There might be something or a combination of things wrong by how we approach SP3 coordinates.
Well, I say wrong 😆 but we keep pushing the limits. So once again, congratulations !

[ChristopherRabotin]

That's encouraging! So at the moment, in that branch, using ANISE is about 10x worse than not using it?

[gwbres]

well it turns out the anise section does not contribute to the kepler solver in the current state, so the previous statement is incorrect.
The last two people's contributions, improved some use cases but made this module a complete mess. Impossible to understand what's happening.

These are the four cases that should be implemented : https://openrtk.readthedocs.io/en/latest/algorithms/EphemeridesClocks.html

[gwbres]

FYI, the initial calculations I implemented come from the Thesis in the bibliography,

Currently, everything that works is based on hepler.orbit_position which is said to be the position in orbital plane.
Any attempts to replace that by Orbit::try_keplerian does not work.
I just added more comments to all these quantities in the code

[ChristopherRabotin]

C'est assez surprenant. Je me demande si les signaux GNSS doivent exister dans un plan orbital spécifique qui n'est pas souvent utilisé dans le reste de l'astro. La rotation de la position dans le plan orbital me fait notamment penser à ça.