skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.43k stars 213 forks source link

Propagation seems to be wrong for TLEs with low eccentricity #235

Closed alepinio closed 5 years ago

alepinio commented 5 years ago

Hi,

I'm stuck with the following. I'm doing many orbit calculations for polar earth satellites and I can't understand why, when the eccentricity falls below some value, the perigee is not where it is supposed to be.

OK case: 1) Take the following TLE

1 99999U          19001.00000000 -.00000000  00000-0  00000-0 0 00001
2 99999 090.0000 270.0000 0100000 000.0000 000.0000 14.00000000000000

2) Propagate the TLE to complete nearly an orbit using, say dt = 1 min , and compute distances for each position to the center of the Earth. As right ascension of ascending node is 0.0 and the argument of perigee is also 0.0, perigee position should be a point on the equator. 3) Find the lowest distance, that's the perigee. Plot it. 4) You will find that the perigee happens for a point very nearly to the equator (OK because right ascension of ascending node is 0.0 and the argument of perigee is also 0.0).

NOT OK case: 1) Take the following TLE (only differs in eccentricity from the former, now is 0010000)

1 99999U          19001.00000000 -.00000000  00000-0  00000-0 0 00001
2 99999 090.0000 270.0000 0010000 000.0000 000.0000 14.00000000000000

2) Repeat step 2 from OK case. 2) Repeat step 3 from OK case. 3) For me, perigee is not a point very nearly to the equator! (not OK, given that right ascension of ascending node is 0.0 and the argument of perigee is also 0.0 then the perigee should a point very near to the equator). The lowest distance to the center of the Earth occurs for me at an angle of approx. 60.3 degrees from the ascending node.

What am I doing wrong?

THX, Alejandro

brandon-rhodes commented 5 years ago

I'm not sure — I've always used the numbers from real satellites, I've not tried setting up artificial situations. The first thing I would do would be to make a rotatable 3D plot of a full orbit. It's been a while since I've tried to set up interactive 3D plots in Python, but a search should bring up pages like this one that might help you get started:

https://stackoverflow.com/questions/33436221/displaying-rotatable-3d-plots-in-ipython-or-ipython-notebook

I always try to get a good plot of a situation first, before asking bigger questions, so that I can see a more full picture of what's going on.

alepinio commented 5 years ago

Quick and dirty. but here you are: https://pastebin.com/D57BHB66

Just run:

python test.py

And you will get a 3d interactive plot of the NOT OK case. You will see something like the following:

selection_002

PD: blue stick is just to easy spot the axis x.

brandon-rhodes commented 5 years ago

Thank you! I'll take a look at it this evening. I wonder if one of those Jupyter interactive sliders would let us vary that TLE parameter and see the orbit redraw into different shapes as we adjust it?

alepinio commented 5 years ago

Thank you very much!

I've never used Jupyter, I can't help you there, sorry.

On the other hand, the plot is 3d...so the script would need some effort to make it iterate different eccentricites, set the "camera" to look from a "cool" angle and took an image of it.

I think that just trying different eccentricities by hand (including the one from the OK case) and taking a look of the plots you will get a good picture of what's happening.

This is head -15 test.py (the TLE is hardcoded there...):

#!/usr/bin/env python2
# https://space.stackexchange.com/a/25959

import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy
import scipy.optimize
from skyfield import api

# set artificial TLE (low eccentricity)
line1, line2 = '''
1 99999U          19001.00000000 -.00000000  00000-0  00000-0 0 00001
2 99999 090.0000 270.0000 0010000 000.0000 000.0000 14.00000000000000
'''.strip().splitlines()

PD: the script does not use scipy.optimize, it's a dead import, but to upload a new pastebin just because of that...

alepinio commented 5 years ago

Something I haven't said yet: my choice of inclination and right ascension is totally arbitrary. I've just selected one that was easy for me to understand. Use yours.

kerel-fs commented 5 years ago

I wonder if one of those Jupyter interactive sliders would let us vary that TLE parameter and see the orbit redraw into different shapes as we adjust it?

Yes. :) I've created this notebook Binder to generate the following animation. You get an interactive version of this if you follow the "launch on binder" link.

output

Maybe this visualization can aid at understanding if the observed change in argument of perigee is an actual issue in the propagation or some other (expected?) effect.

brandon-rhodes commented 5 years ago

Thank you for the visualization, @kerel-fs! I've been too busy this spring to try tackling this issue, and I suspect that SGP4 simply does not do well with non-realistic orbits and maybe was not designed to; but I'm happy to see conversation continue.

MichaelRThompson commented 5 years ago

Confirmed the above results with the my own SGP4 propagator. It's not something inherently wrong with the propagator.

SGP4 takes into account multiple perturbations that are messing with your argument of perigee. As your orbit gets closer to circular (a lower eccentricity value), argument of perigee becomes more and more sensitive given that it's closer to a circular orbit, where argument of perigee is pretty much undefined.

For your two test cases:

You're approaching your problem with the assumption that your argument of perigee is going to stay relatively constant, but in reality it is varying by quite a bit. Your two test cases are at about 800 km altitude, and the perturbations due to things like J2 are still quite strong there.

I created a 3rd test case using the TLE:

1 99999U          19001.00000000 -.00000000  00000-0  00000-0 0 00001
2 99999 090.0000 270.0000 0010000 000.0000 000.0000 5.00000000000000

This is a much higher orbit, and the argument of perigee is bounded between ~5 and 45 degrees.

As you can see, your issue isn't really a problem with the propagator, it's an issue with perturbations. Especially at low eccentricities, your argument of perigee can become really sensitive to perturbations, and your periapsis can move around a good bit.

brandon-rhodes commented 5 years ago

Thank you for the analysis, @MichaelRThompson! Can anyone suggest an improvement to the documentation that might cover this case? For now I'll close the issue, but if anyone wants to submit some updated docs, please feel free to do so; otherwise, hopefully Google will lead folks to this issue who have a similar question.

MichaelRThompson commented 5 years ago

I'm not sure that this case specifically needs to be addressed, but looking over the documentation, a quick introduction to the SGP4 propagator could be helpful, probably on the Earth Satellites page here. Just an overview of what kind of perturbations are included and what that can do to some of your orbital elements.

I think it could fit it nicely with the limitations section up at the top. Let me know if you'd like any help in writing something up.