Fix near parabolic regimes

[astrojuanlu]

Regions of (nu, ecc) from Farnocchia at al., 2013:

Previously our logic was wrong: we were incorrectly using the near parabolic solution where we should have been using the strong elliptic solution, and therefore it failed "in the vicinity of nu = pi, 0.985 < ecc < 1.0 (for delta=1e-2)" as I correctly identified two years ago.

• Fix nu_to_M (which causes orbit.M to hang)
• Fix M_to_nu (which causes the mean_motion propagator to hang)

Fix #475 (finally!), fix #907.

[astrojuanlu]

Asked one question about notation in Math.StackExchange - it turns out that T^1 is the 1-torus, the circle. Cool!

[codecov]

Codecov Report

Merging #908 into master will decrease coverage by 0.59%.
The diff coverage is 76.41%.

@@            Coverage Diff             @@
##           master     #908      +/-   ##
==========================================
- Coverage   89.93%   89.34%   -0.60%     
==========================================
  Files          69       69              
  Lines        3608     3668      +60     
  Branches      301      316      +15     
==========================================
+ Hits         3245     3277      +32     
- Misses        286      308      +22     
- Partials       77       83       +6     

Impacted Files	Coverage Δ
src/poliastro/twobody/orbit.py	87.02% <72.00%> (-1.23%)	⬇️
src/poliastro/core/propagation/farnocchia.py	79.02% <77.63%> (-10.88%)	⬇️
src/poliastro/core/angles.py	95.29% <80.00%> (+0.17%)	⬆️
src/poliastro/twobody/states.py	80.73% <0.00%> (-0.92%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 2a01615...7f5dc6b. Read the comment docs.

[astrojuanlu]

Well, this was more difficult than expected. Getting errors in some old tests that compare the propagator with Cowell's method, still not sure what's happening. Will review tomorrow.

[astrojuanlu]

Well, this turned out to be more difficult than I thought.

Several findings while working on this:

• Fixing the near parabolic region computation in nu_to_M and M_to_nu fixes some cases, but breaks others. The functions does not seem to be smooth at the boundaries, so small numerical errors can have a big effect and produce totally different results.
• We have a huge mess in the Orbit.sample() method and it's urgent that we refactor it to avoid several useless conversions between classical and cartesian elements.
• There is a very tight coupling between how the Farnocchia method separates the "near parabolic" region, and all the code in poliastro.core.angles. Which either means:
  • There shouldn't be such coupling, and all the near-parabolic code should be used only by poliastro.core.propagation.farnocchia.
  • All the propagators need such distinctions to work properly for near-parabolic orbits, and therefore the near-parabolic boundary will trascend the farnocchia propagator.

I still don't understand why nu_to_M and M_to_nu are not smooth. There's excellent mathematical background about them in the paper, but it requires deeper inspection.

We have a path to fix this, but it will take a long time.

[astrojuanlu]

Paper	Numerical simulation

Looking great so far! Now on to the next step 💪

[astrojuanlu]

All the propagators in poliastro.core.propagation use some sort of angle conversion at some point, and in some algorithms branch depending on the eccentricity.

A very obvious case can be seen in Mikkola's method:

poliastro/src/poliastro/core/propagation.py

Lines 327 to 335 in 643356d

	if ecc < 1.0:
	nu = E_to_nu(E, ecc)
	else:
	if ecc == 1.0:
	# Parabolic
	nu = D_to_nu(E)
	else:
	# Hyperbolic
	nu = F_to_nu(E, ecc)

I feel that if we mix the treatment that Farnocchia does of the near-parabolic region with the other propagators we can get mixed results. Doing so probably requires a careful review of all the algorithms and implementations. The last two years we have devoted a lot of time to improve our propagators, and I have to say both @jorgepiloto and @nikita-astronaut did an amazing job.

Therefore, I will try not to affect how other propagators work while fixing the issues with the Farnocchia one.

And to finally justify all this work, let's take a look at this table again:

poliastro/src/poliastro/twobody/propagation.py

Lines 4 to 22 in 643356d

	+-------------+------------+-----------------+-----------------+
	| Propagator | Elliptical | Parabolic | Hyperbolic |
	+-------------+------------+-----------------+-----------------+
	| farnocchia | ✓ | ✓ | ✓ |
	+-------------+------------+-----------------+-----------------+
	| vallado | ✓ | ✓ | ✓ |
	+-------------+------------+-----------------+-----------------+
	| mikkola | ✓ | NOT IMPLEMENTED | ✓ |
	+-------------+------------+-----------------+-----------------+
	| markley | ✓ | x | x |
	+-------------+------------+-----------------+-----------------+
	| pimienta | ✓ | ✓ | NOT IMPLEMENTED |
	+-------------+------------+-----------------+-----------------+
	| gooding | ✓ | NOT IMPLEMENTED | NOT IMPLEMENTED |
	+-------------+------------+-----------------+-----------------+
	| danby | ✓ | x | ✓ |
	+-------------+------------+-----------------+-----------------+
	| cowell | ✓ | ✓ | ✓ |
	+-------------+------------+-----------------+-----------------+

Leaving Cowell's method aside, Farnocchia's and Vallado's are the only propagators we have at the moment that are supposed to work with all orbits. Vallado's method, being based on the universal variables approach, shares its pros and cons. From Farnocchia et al.:

Even though @nikita-astronaut fixed the convergence problems of Vallado's method in #362, it still requires a large number of Newton iterations (in line with the comment above) and the relative tolerance is bounded by numerical errors.

In conclusion:

• It's worth it to try to fix this propagator, and I will do so in this PR.
• The angles functions should be work as before because they are used in other propagators.
• We will keep Farnocchia's method as the one preferred for sampling, and we will do some refactoring to avoid unnecessary element conversions and accelerate it considerably.

[astrojuanlu]

While working on this PR I came to an important realization:

Even if most books and papers talk about the mean anomaly as a quantity that can be defined for elliptic, parabolic and hyperbolic orbits, more rigorous texts like Battin are very careful not to abuse the term for eccentricities equal or larger than one.

For example, the "mean anomaly" for hyperbolic motion is referred to as "the quantity N":

The same can be said about the "parabolic anomaly" and the "hyperbolic anomaly", which are never called as such: the former is mentioned tangentially suggested when saying that "tan f/2 is the key variable in Barker's equation", and the latter is called "the quantity H", or the "hyperbolic analog to the eccentric anomaly".

My theory is that this is done because these quantities do not represent angles directly. However, for practical reasons this fact is often overlooked.

[astrojuanlu]

My previous comment has important implications for poliastro. In the same way that Orbit objects do not have an E property (for the eccentric anomaly), they shouldn't have an M property for the mean anomaly either, and it should be treated at a lower level as a useful mathematical artifact, but nothing more.

[astrojuanlu]

After the initial guess change, the following tests are still failing on CI:

• test_long_propagation_preserves_orbit_elements[farnocchia] (macOS)
• test_long_propagations_vallado_agrees_farnocchia (macOS)

The Windows build passed though!

[helgee]

Those pass locally with NUMBA_DISABLE_JIT=1 on my machine. But fail while jitted.

[astrojuanlu]

Can't reproduce these locally, either with numba 0.48 or 0.49, with or without jitting.

[noc0lour]

Those pass locally with NUMBA_DISABLE_JIT=1 on my machine. But fail while jitted.

It would be interesting to know the processor type/supported SIMD instruction sets, if numba is optimizing w.r.t. those that might be a difference...

[helgee]

It would be interesting to know the processor type/supported SIMD instruction sets, if numba is optimizing w.r.t. those that might be a difference...

Intel(R) Core(TM) m7-6Y75 CPU @ 1.20GHz

https://ark.intel.com/content/www/us/en/ark/products/88199/intel-core-m7-6y75-processor-4m-cache-up-to-3-10-ghz.html

[astrojuanlu]

Found another falsifying case locally: nu_from_delta_t(3393554832.0050254, 0.9676000001129921, 132712442099.00002, 87574593.50447243)

[astrojuanlu]

All tests passed 😍

Thanks everyone who chimed in! I am now merging this once and for all.

[noc0lour]

@helgee are your failing tests (with JIT) covered with the fixes in the latest commits? If not, it might be worth to open a separate issue on these.

[helgee]

Works fine now 👍

[mkbrewer]

@astrojuanlu Just a followup on your comment regarding issue #10230 in the astropy repository. I'm posting here as it seems like a more appropriate place to do so. Your comment was a bit concerning as I have a program that is in use by the control systems at several observatories to track comets in real time given a set of orbital elements. The algorithm that I use is attributed to Danby from the 2nd edition of his book Fundamentals of Celestial Mechanics 1988 (pages 149-154). I don't have a copy of this book, but it uses an accelerated Newton's method very similar to the one I see in poliastro/core/propigation.py. The only difference that I see is that:

  fpp = s / 2.0
  fppp = c / 6.0

I have never had any problems with this program. From what I see here, the issues that you were seeing only occur when the true anomaly is close to 180 degrees. Is that correct? If that is the case, then I'm not going to worry much as we are only concerned with observing comets as they make their closest approach to Earth.

[mkbrewer]

Oh, never mind on the difference in the algorithms. I see that yours has the factors of 1/2 and 1/6 absorbed in to deltas whereas mine doesn't.

[astrojuanlu]

Thanks for chiming in @mkbrewer! (xref astropy/astropy#10230)

The conversation in this particular pull request is a mess. The important thing is the end result, based on the Farnocchia et al. paper. Let me know if you have any comments about it.

[mkbrewer]

No comments, just my question above.

[mkbrewer]

@astrojuanlu An answer to this question would be appreciated when you have the time.

From what I see here, the issues that you were seeing only occur when the true anomaly is close to 180 degrees. Is that correct?

[astrojuanlu]

Hello @mkbrewer, I'm so sorry for the huge delay - I was in a bad moment when you asked, and then when I recovered I have been postponing it. Trying to get back on track.

From what I see here, the issues that you were seeing only occur when the true anomaly is close to 180 degrees. Is that correct?

Yes, that's correct. The issues that we were seeing in this algorithm were for the case where the true anomaly is close to 180 degrees. If you think about it, that condition is impossible to happen in true parabolic motion, and for near-parabolic motion, the closer we are to 180 the more evidence we have that the orbit is, well, elliptic. That's why in the chart by Farnocchia et al. the "near parabolic" region shrinks towards 180 degrees:

Fortunately, we managed to fix our implementation (it contained several conceptual errors) and now our farnocchia propagator can handle orbits with eccentricity arbitrarily close to 1.0, and do a smooth transition in the frontiers of the regions.

Hope this answers your question, please let me know otherwise.