shashwatak / satellite-js

Modular set of functions for SGP4 and SDP4 propagation of TLEs.
MIT License
902 stars 145 forks source link

Two Errors in Existing Implementation #107

Closed TSKelso closed 1 year ago

TSKelso commented 2 years ago

In performing a full-catalog (23,781 TLEs) comparison of satellite.js to the official results from the US Space Force official results from SGP4Prop, I discovered two issues with the current implementation.

  1. The current version incorrectly chooses the WGS84 set of geodetic parameters. While it is a more current geodetic system than WGS72, USSF still uses WGS72. It is important to use the SAME version of SGP4 used to produce the TLE (GP) data when propagating. For a test case using SSC 00005, using WGS84 values produced an 11.5-m difference. Making the change below resulted in a difference of zero (to 6 decimal places).

    / WRONG VALUES (WGS84) var mu = 398600.5; // in km3 / s2 var earthRadius = 6378.137; // in km var j2 = 0.00108262998905; var j3 = -0.00000253215306; var j4 = -0.00000161098761; / // CORRECT VALUES (WGS72) var mu = 398600.8; // in km3 / s2 var earthRadius = 6378.135; // in km var j2 = 0.001082616; var j3 = -0.00000253881; var j4 = -0.00000165597;

  2. Even after making this change, there were 377 cases out of 23,781 where the differences were non-zero. Inspection of those cases revealed they all involved eccentric orbits (e > 0.5) that were nearly semi-synchronous (mean motion ~2 revs/day). Comparing the C# version of the code in the section labeled "near - half-day resonance terms" to that in satellite.js revealed a missing set of parentheses on the 6th line of that section. The missing parentheses are shown with the # symbols in the corrected equation below:

xnddt = d2201 Math.cos(x2omi + xli - g22) + d2211 Math.cos(xli - g22) + d3210 Math.cos(xomi + xli - g32) + d3222 Math.cos(-xomi + xli - g32) + d5220 Math.cos(xomi + xli - g52) + d5232 Math.cos(-xomi + xli - g52) + 2.0 #(d4410 Math.cos(x2omi + x2li - g44) + d4422 Math.cos(x2li - g44) + d5421 Math.cos(xomi + x2li - g54) + d5433 * Math.cos(-xomi + x2li - g54))#;

Once this correction is made, propagation of all 23,781 TLEs with satellite.js agrees with the USSF results from the official version of SGP4Prop (Version 8.3) to 6 decimal places (i.e., less than 1 mm).

thkruz commented 2 years ago

I saw in Revisiting Spacetrack Report 3 that the assumption was made that (then) AFSPC would likely stick to WGS-72 instead of upgrading because of interoperability issues on older systems that aren't updated. That seemed like a fair assumption and my fork that I hope will one day replace this library defaults to WGS-72, but since the creator of this repo has moved on, I have been apt to make serious changes to the code.

I have zero doubt in my mind that you are correct about the section on half-day resonance terms. Downloading SGP4Prop to generate a few test cases to ensure this doesn't get broke again.

@ezze now that I got my TypeScript version working I have completely stopped using this package - are you still around and able to merge a PR if I make it or is this one more reason to deprecate this package?

ezze commented 2 years ago

@thkruz Hi, mate! I am still available, and it would be nice if we could upgrade the package because many people still use it. Looking on how repository's stars are growing we can be sure that our work is still in demand.

I see that changes are not too complicated but we must make them carefully not breaking the tests. We have not too many tests, and I'm not sure whether at least one exists to cover the case described above. So if it's possible please also check this.

thkruz commented 2 years ago

I'll take another stab after work. I found satellite 692 that will trigger that block of code:

"tleLine1": "1 00692U 63039C   22222.54784188 -.00000483  00000+0  00000+0 0  9997",
"tleLine2": "2 00692  10.2323 330.8569 5575296 261.9869   0.5612  1.99110435 47456",

I can see in dspace.js where xnddt is different with/without the second modification @TSKelso recommended, but it isn't having any impact on the position/velocity output.

As for the WGS change - that will break every single position/velocity check in the tests because they are currently incorrect. Ideally we would test against the "official" values in SGP4Prop's test. I'll upload those and then it is just a matter of putting them in the json format.

ezze commented 2 years ago

Ok, then we must double check updated fixtures for the tests.

TSKelso commented 2 years ago

You should not be testing with cases from Spacetrack Report Number 3. That was initially published in 1980 and many errors have been detected and corrected in that code. Those tests only cover a couple of cases, anyway.

You should be using the USSF official SGP4, which is available on Space Track at (you need a Space Track account):

https://www.space-track.org/documentation#/sgp4 <https://www.space-track.org/documentation#/sgp4>

I used the latest version (Version 8.3) to produce results for all 23,781 objects we had TLE data for on Jul 20 for a period of ±1440 minutes from TLE epoch at 1-minute intervals for a test case. I then set up a similar test with satellite.js. All positions were reported to 6 decimal places or to the 1-mm level.

Doing that with the current version of satellite.js showed errors at the meter level for all satellites, meaning the error was something that affected all orbits. I started at the top of the code where I noticed that the Earth radius and Earth gravitational parameter were not the WGS72 values. Many people assume that since WGS84 is an updated version of WGS72, that it should be used instead. However, USSF has never made that change. Once I replaced the WGS84 parameters with WGS72 ones, all but 377 cases had their differences go to zero (0.000000) km.

In all 377 cases, those were eccentric semi-synchronous orbits—not one of your test cases. That allowed me to find the appropriate section of the code and compare it to the code from our paper, which is available on CelesTrak. Doing a text comparison clearly showed the missing parentheses. Correcting that oversight resulted in all 23,781 cases matching to 0.000000 km.

BTW, I had also done a full comparison of the C# code from our paper to the USSF SGP4 and found full agreement for all 23,781 cases (so perfect I thought I must have performed the test incorrectly). So, that clearly verifies that USSF is still using WGS72.

It took me about two weeks of focused testing to implement both SGP4Lib (our code), SGP4Prop (USSF version), and satellite.js. I wanted to be sure not just to pick a couple of random test cases that happened to work and then move on. I would likely not have found the missing parentheses if I hadn’t tested the full catalog, since 377 of 23,781 is only 1.6% of all cases.

I am sure you can verify these cases independently, but you might expect to spend a similar amount of time in validating and testing each piece of code. Just let me know if you have any questions. - TS

Dr. T.S. Kelso CelesTrak, https://celestrak.org E-Mail: @.***

See how you can help support CelesTrak’s non-profit educational mission by donating at https://celestrak.org today!

On 2022 Aug 09, at 01:24, Dmitriy Pushkov @.***> wrote:

Ok, then we must double check updated fixtures for the tests.

— Reply to this email directly, view it on GitHub https://github.com/shashwatak/satellite-js/issues/107#issuecomment-1209254796, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALCCCZTJVKTMWBBEQNI42MDVYI5XPANCNFSM55XXKJEA. You are receiving this because you were mentioned.

thkruz commented 2 years ago

@TSKelso Copy. The explanation is immensely helpful.

Downloaded SGP4Prop and generated position/velocity data for ~26,000 satellites at 1 minute intervals for 1440 minutes. Then wrote a node script to parse the output from SGP4Prop and format it in JSON. Jest kept barfing when trying to test every minute even when split into groups of 1000 satellites. I reduced the testing to checking at minutes 0, 360, 720, 1080, and 1440 and I am not getting 0.000000km accuracy after making the two changes (or before for that matter).

When looking only at epoch (+0 minutes) 92% of the satellites are within 0.000005km

All but the following two are within 0.00005km:

1  4867U 70025EX  22221.79749997  .00000482  00000-0  33778-3 0  9994
2  4867  99.6797  69.1653 0136587 331.8508  27.5317 13.99383662658130
1 46657U 18079CK  22221.26010476 -.00000034  00000-0  00000+0 0  9995
2 46657   8.9156 167.9825 5390259  63.9254 341.0492  1.69839457 11387

And they are off by 0.05 meters and 0.07 meters respectively. For my use case that is more than close enough. Unless you are calculating probability of collision or trying to dock, not sure a fraction of a meter is going to matter, but I am very intrigued that I am not getting the same accuracy if we are both using satellite.js and SGP4Prop 8.3. I am using Node 17.2.0 on Windows 11 for satellite.js and used the Python wrapper in SGP4Prop to get my truth data. I will try it against ootk tomorrow since that one more closely mirrors the C# code your paper.

Now when looking at the 360, 720, 1080, and 1440 minute markers I am getting everything within 0.5 meters at 360 and 720 and within 5 meters at 1080 and 1440 except:

1 48357U 21038E   22221.42821302  .00004990  00000-0  35347-3 0  9999
2 48357  53.0538 114.6689 0001852  76.1103 284.0092 15.06419300 70430
1 32984U 79101K   22219.47995179 -.00000019  00000-0  00000-0 0  9992
2 32984   8.5154 223.4941 4863823 197.4941 284.0078  1.82638276 63927
1 82117U          18267.39743894 +.00001817 +00000-0 +74703-3 0  9991
2 82117 069.9048 280.9712 0017481 316.1290 043.8443 14.28485864747744

Which are often off by 1-10km especially 82117. Additionally this one:

1 50298U 82092ZP  22221.74953100  .00112912  00000-0  18977-2 0  9993
2 50298  82.8796 227.5408 0037647 234.2150 125.5602 15.50541526 39218

Fails to calculate position/velocity at 1440 minutes. I will see if I can find what is causing the error.

The other 130300+ tests are within tight margins. So I will look at those edge cases tomorrow, but the overwhelming majority are showing up really close to SGP4Prop, just not as close as @TSKelso was getting them.

@ezze the one giant issues I am running into is that the output of SGP4Prop comes with a disclaimer:

*  THIS SOFTWARE CONTAINS TECHNICAL DATA WHOSE EXPORT IS *
*  RESTRICTED BY THE ARMS EXPORT CONTROL ACT (TITLE 22,  *
*  USC, SEC 2751 ) OR EXECUTIVE ORDER 12470. VIOLATORS   *
*  OF EXPORT LAWS ARE SUBJECT TO SEVERE CRIMINAL         *
*  PENALTIES.

So there is no way I am uploading the test cases (not to mention they are currently sitting at 5 Gb of json). How would you feel about me using snapshots once we are happy with the results? That way you will have immediate feedback if any changes impact the accuracy (over 100,000 test cases) without the need to know the exact values USSF generated.

TSKelso commented 2 years ago

Theodore,

That’s odd. I had no problem generating that data. I used a full dataset in TLE format from CelesTrak (from Jul 20, so all tests used the same data). In SGP4Lib (the C# code from our 2006 paper), I looped through each TLE, generated data ±1440 minutes from TLE epoch at 1-minute intervals, and output that to a file. That prevents using a lot of internal memory, since it is just propagating one satellite at a time. I did the same thing for SGP4Prop. I then did the same thing in satellite.js, except that I added the comparison step.

As indicated, no data agreed using WGS84 coefficients instead of WGS72. Changing that got 98.4% agreement. Being able to search for any non-zero results (I added “***” as a search key for those cases) allowed me to review cases that failed and see the pattern for the type of orbit, which led to the section of the code. A quick comparison to the SGP4Lib C# code showed the missing parentheses and eliminated all discrepancies. There is no doubt that those parentheses are missing from satellite.js.

To avoid the issue with the Arms Export Control Act, you can simply validate against the SGP4Lib results, which do not have that restriction. Of course, if you are not sure that version is working correctly, you can compare it to SGP4Prop, as I did.

The JavaScript in satellite.js should match SGP4Lib, not a Python version of SGP4Lib that may have the same issues that I found in satellite.js.

For my testing, the SGP4Lib DLL was built using Visual Studio 2022 (Version 17.2) as 64-bit. Both the SGP4Lib and SGP4Prop programs to produce propagated output were run in that environment, too. Everything was done on Windows Server 2019. The satellite.js test was performed on the same Server 2019 platform running IIS.

Just let me know if you have any questions. - TS

Dr. T.S. Kelso @.***

On 2022 Aug 09, at 17:50, Theodore Kruczek @.***> wrote:

@TSKelso https://github.com/TSKelso Copy. The explanation is immensely helpful.

Downloaded SGP4Prop and generated position/velocity data for ~26,000 satellites at 1 minute intervals for 1440 minutes. Then wrote a node script to parse the output from SGP4Prop and format it in JSON. Jest kept barfing when trying to test every minute even when split into groups of 1000 satellites. I reduced the testing to checking at minutes 0, 360, 720, 1080, and 1440 and I am not getting 0.000000km accuracy after making the two changes (or before for that matter).

When looking only at epoch (+0 minutes) 92% of the satellites are within 0.000005km

All but the following two are within 0.00005km:

1 4867U 70025EX 22221.79749997 .00000482 00000-0 33778-3 0 9994 2 4867 99.6797 69.1653 0136587 331.8508 27.5317 13.99383662658130 1 46657U 18079CK 22221.26010476 -.00000034 00000-0 00000+0 0 9995 2 46657 8.9156 167.9825 5390259 63.9254 341.0492 1.69839457 11387 And they are off by 0.05 meters and 0.07 meters respectively. For my use case that is more than close enough. Unless you are calculating probability of collision or trying to dock, not sure a fraction of a meter is going to matter, but I am very intrigued that I am not getting the same accuracy if we are both using satellite.js and SGP4Prop 8.3. I am using Node 17.2.0 on Windows 11 for satellite.js and used the Python wrapper in SGP4Prop to get my truth data. I will try it against ootk tomorrow since that one more closely mirrors the C# code your paper.

Now when looking at the 360, 720, 1080, and 1440 minute markers I am getting everything within 0.5 meters at 360 and 720 and within 5 meters at 1080 and 1440 except:

1 48357U 21038E 22221.42821302 .00004990 00000-0 35347-3 0 9999 2 48357 53.0538 114.6689 0001852 76.1103 284.0092 15.06419300 70430 1 32984U 79101K 22219.47995179 -.00000019 00000-0 00000-0 0 9992 2 32984 8.5154 223.4941 4863823 197.4941 284.0078 1.82638276 63927 1 82117U 18267.39743894 +.00001817 +00000-0 +74703-3 0 9991 2 82117 069.9048 280.9712 0017481 316.1290 043.8443 14.28485864747744 Which are often off by 1-10km especially 82117. Additionally this one:

1 50298U 82092ZP 22221.74953100 .00112912 00000-0 18977-2 0 9993 2 50298 82.8796 227.5408 0037647 234.2150 125.5602 15.50541526 39218 Fails to calculate position/velocity at 1440 minutes. I will see if I can find what is causing the error.

The other 130300+ tests are within tight margins. So I will look at those edge cases tomorrow, but the overwhelming majority are showing up really close to SGP4Prop, just not as close as @TSKelso https://github.com/TSKelso was getting them.

@ezze https://github.com/ezze the one giant issues I am running into is that the output of SGP4Prop comes with a disclaimer:

  • THIS SOFTWARE CONTAINS TECHNICAL DATA WHOSE EXPORT IS *
  • RESTRICTED BY THE ARMS EXPORT CONTROL ACT (TITLE 22, *
  • USC, SEC 2751 ) OR EXECUTIVE ORDER 12470. VIOLATORS *
  • OF EXPORT LAWS ARE SUBJECT TO SEVERE CRIMINAL *
  • PENALTIES. So there is no way I am uploading the test cases (not to mention they are currently sitting at 5 Gb of json). How would you feel about me using snapshots once we are happy with the results? That way you will have immediate feedback if any changes impact the accuracy (over 100,000 test cases) without the need to know the exact values USSF generated.

— Reply to this email directly, view it on GitHub https://github.com/shashwatak/satellite-js/issues/107#issuecomment-1210126138, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALCCCZV2HVFHYRS6EGUKN2LVYMRJFANCNFSM55XXKJEA. You are receiving this because you were mentioned.

ezze commented 1 year ago

Fixed in 5.0.0.