Wrong predictions for Molniya orbits

[emedvedev]

I'm trying to predict the position of a satellite with a Molniya orbit, Molniya 3-47:

1 23642U 95042A   18058.62866679  .00000210  00000-0 -75405-2 0  9996
2 23642  62.9606 257.1441 7218342 288.9618  10.0525  2.01983768165381

The prediction turns out to be completely incorrect. The minimal example below is taken from the readme almost line-by-line:

var tleLine1 = '1 23642U 95042A   18058.62866679  .00000210  00000-0 -75405-2 0  9996',
    tleLine2 = '2 23642  62.9606 257.1441 7218342 288.9618  10.0525  2.01983768165381';
var satrec = twoline2satrec(tleLine1, tleLine2);
var positionAndVelocity = propagate(satrec, new Date());
var positionEci = positionAndVelocity.position;
var gmst = gstime(new Date());
var positionGd = eciToGeodetic(positionEci, gmst);
console.log('Molniya', positionGd.latitude * (180 / Math.PI), positionGd.longitude * (180 / Math.PI));

It gives me 39.8491516985329 -157.00043723551485, while the correct coordinates are 38.64 72.56. At other times it's even more off. Here's the plotting of the orbit on Google Maps:

In reality it's supposed to look like this (see n2yo):

Predictions for any other orbit (including highly elliptical orbits, like the Tundra orbit of QZSS satellites) seem to be correct.

[emedvedev]

A correct result from pyephem:

>>> import ephem
>>> line1 = 'Molniya 3-47'
>>> line2 = '1 23642U 95042A   18058.62866679  .00000210  00000-0 -75405-2 0  9996'
>>> line3 = '2 23642  62.9606 257.1441 7218342 288.9618  10.0525  2.01983768165381'
>>> molniya = ephem.readtle(line1, line2, line3)
>>> molniya.compute()
>>> molniya.sublong
76:15:28.8
>>> molniya.sublat
45:42:17.3

[ezze]

Hi, @emedvedev. Thank you for the report. I confirm that something is broken in 2.0.0. I get different results with 1.3.0 and they are close to expected. We need some time to explore what's going wrong.

[ezze]

After some investigation I stated that results' difference between versions 1.3.0 and 2.0.0 is caused by this pull request that fixes #26. To be exact, it's a fix in dsinit function:

   // -------------------- deep space initialization ------------
   irez = 0;
-  if (nm > 0.0034906585 < 0.0052359877) { // eslint-disable-line
+  if ((nm < 0.0052359877) && (nm > 0.0034906585)) {
     irez = 1;
   }
-  if (nm >= 8.26e-3 <= 9.24e-3 && em >= 0.5) { // eslint-disable-line
+  if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5)) {
     irez = 2;
   }

So in 1.3.0 this code snippet results in irez = 0, in 2.0.0 — irez = 2. If you set irez to 0 then you will get the expected calculations.

Unfortunatelly, I have no opportunity to make similar calculations with original SGP4 python library right now but, probably, @tikhonovits, @dangodev, @shashwatak can look at this. What do you think, guys?

[drwpow]

For some context, the same part of the SGP4 code, compared across libraries:

satellite-js

// -------------------- deep space initialization ------------
irez = 0;
if ((nm < 0.0052359877) && (nm > 0.0034906585)) {
  irez = 1;
}
if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5)) {
  irez = 2;
}

Original SPG4 Code (C++)

/* -------------------- deep space initialization ------------ */
irez = 0;
if ((nm < 0.0052359877) && (nm > 0.0034906585))
    irez = 1;
if ((nm >= 8.26e-3) && (nm <= 9.24e-3) && (em >= 0.5))
    irez = 2

Python sgp4 (1.4)

#  -------------------- deep space initialization ------------
irez = 0;
if 0.0034906585 < nm < 0.0052359877:
    irez = 1;
if 8.26e-3 <= nm <= 9.24e-3 and em >= 0.5:
    irez = 2;

As @ezze mentioned, the 2.0.0 bump for satellite-js was restoring the values from the original C++ that had been changed in the Python port. After the change, I saw much improved accuracy in ISS location-tracking for my usecase, as did @tikhonovits in the original issue that led to the change.

As to why they were changed in Python’s version—which there’s probably a good reason—I can’t really say.

[nhamer]

nb, the orbit for QZS-1 is irez=1 (1 day resonance, rather than 2=half-day)

1 37158U 10045A 18103.52828076 -.00000124 00000-0 00000+0 0 9992 2 37158 41.0319 154.9738 0753103 270.0506 104.4433 1.00281986 27777

[nhamer]

Almost certain that the issue is https://github.com/shashwatak/satellite-js/blob/226d8d5520c0cf03a591a6f0bd4ddfd8a1de0965/src/propagation/dspace.js#L248

The equivalent C code (sgpunit.cpp@1114) is: mm = xl - 2.0 * nodem + 2.0 * theta;

(slightly different runs over a few seconds, so gmstime is not identical):

Old satellite.js (irez=0): Molniya 62.97302163853729 -152.31526401326246 Current satellite.js (irez=2): Molniya 30.900731334965187 -76.42315864995993 Correcting mm irez=2) gives: Molniya 62.864247526501345 -155.86674453746895

[ezze]

@nhamer, nice catch! Thanks a lot for this.

Following Airbnb code style forced me to add many parentheses to satisfy no-mixed-operators rule. As we see now, it doesn't give any benefits but adds headaches like this issue. Probably, it would be better to disable this rule in .eslintrc.json as it was done for function-paren-newline later. It was my mistake, sorry for this, guys.

A fix is made in https://github.com/shashwatak/satellite-js/commit/d398bd5fb880a48708c3ef6be712a595360c12b3. You can update the library to 2.0.2. @emedvedev, please check whether everything works fine for you. @dangodev, hope nothing is broken for you too.

All tests are passed but we don't have too much tests to cover everything. If someone of you, guys, can prepare some test fixtures for Molniya 3-47 with C++ or Python version of the library, you are welcome. :)

[emedvedev]

@ezze thanks! Seems to be working now, and yes, default code styles can be a pain sometimes. :)